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We review the present status of our research and understanding regarding the dynamics and the 
statistical properties of earthquakes, mainly from a statistical physical viewpoint. Emphasis is 
put both on the physics of friction and fracture, which provides a "microscopic" basis of our 
understanding of an earthquake instability, and on the statistical physical modelling of earth- 
quakes, which provides "macroscopic" aspects of such phenomena. Recent numerical results on 
several representative models are reviewed, with attention to both their "critical" and "character- 
istic" properties. We highlight some of relevant notions and related issues, including the origin of 
power-laws often observed in statistical properties of earthquakes, apparently contrasting features 
of characteristic earthquakes or asperities, the nature of precursory phenomena and nucleation 
processes, the origin of slow earthquakes, etc. 
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I. INTRODUCTION 

Earthquakes are large scale mechanical failure phenom- 
ena, which have still defied our complete understand- 
ing. In this century, we already experienced two gi- 
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gantic earthquakes: 2004 Sumatra- Andaman earthquake 
(M9.1) and 2011 East Japan earthquake (M9.0). Given 
the disastrous nature of the phenomena, the understand- 
ing and forecasting of earthquakes have remained to be 
the most important issue in physics and geoscience (Carl- 
son, Langer and Shaw, 1996; Rundle, Turcotte and Klein, 
2000; Scholz, 2002; Rundle et al, 2003; Bhattacharyya 
and Chakrabarti, 2006; Ben-Zion, 2008; Burridge, 2006; 
De Rubies et al, 2006; Kanamori, 2009; Daub and Carl- 
son, 2010). Although there is some recent progress in 
our understanding of the basic physics of fracture and 
friction, it is still at a primitive stage (Marone, 1998; 
Scholz, 1998, 2002; Dieterich, 2009; Tulhs, 2009; Daub 
and Carlson, 2010). Furthermore, our lack of a proper 
understanding of the dynamics of earthquakes poses an 
outstanding challenge to both physicists and seismolo- 
gists. 

While earthquakes are obviously complex phenomena, 
certain empirical laws have been known concerning their 
statistical properties, e.g., the Gutenberg- Richter (GR) 
law for the magnitude distribution of earthquakes, and 
the Omori law for the time evolution of the frequency of 
aftershocks (Scholz, 2002; Rundle, 2003; Turcotte, 2009). 
The GR law states that the frequency of earthquakes of 
its energy (seismic moment) E decays with E obeying 
a power-law, i.e., oc £'~(i+^) = £;~(i+t'') where B and 
b = |-B are appropriate exponents, whereas the Omori 
law states that the frequency of aftershocks decays with 
the time elapsed after the mainshock obeying a power- 
law. These laws, both of which are power-laws possess- 
ing a scale- invariance, are basically of statistical nature, 
becoming evident only after examining large number of 
events. Although it is extremely difficult to give a defini- 
tive prediction for each individual earthquake event, clear 
regularity often shows up when one measures its statis- 
tical aspect for an ensemble of many earthquake events. 
This observation motivates statistical physical study of 
earthquakes due to the following two reasons: First, a 
law appearing after averaging over many events is ex- 
actly the subject of statistical physics. Second, a power 
law or a scale invariance has been a central subject of 
statistical physics over years in the context of critical 
phenomena. Indeed, Bak and collaborators proposed 
the concept of "self-organized criticality (SOC)" (Bak, 
Tang and Wiesenfeld, 1987). According to this view, 
the Earth's crust is always in the critical state which 
is self-generated dynamically (Turcotte 1997; Hergarten, 
2002; Turcotte, 2009; Pradhan, Hansen and Chakrabarti, 
2010). One expects that such an SOC idea might possi- 
bly give an explanation of the scale-invariant power-law 
behaviors frequently observed in earthquakes, including 
the GR law and the Omori law. However, one should 
also bear in mind that real earthquakes often exhibit ap- 
parently opposite features, i.e., the features represented 
by "characteristic earthquakes" where an earthquake is 
regarded to possess its characteristic energy or time scale 
(Scholz, 2002; Turcotte, 2009). 

Earthquakes also possess strong relevance to material 



science. It is now established that earthquakes could 
be regarded as a stick-slip frictional instability of a pre- 
existing fault, and statistical properties of earthquakes 
are governed by the physical law of rock friction (Marone, 
1998; Scholz, 1998; 2002, Dieterich, 2009; Tullis, 2009). 
The physical law describing rock friction or fracture is 
often called "constitutive law". As most of the major 
earthquakes are caused by rubbing of faults, such fric- 
tion laws give the "microscopic" basis in analyzing the 
dynamics of earthquakes. One might naturally ask: How 
statistical properties of earthquakes depend on the mate- 
rial properties characterizing earthquake faults, e.g., the 
elastic properties of the crust or the frictional properties 
of the fault, etc. Answering such questions would give 
us valuable information in understanding the nature of 
earthquakes. 

In spite of some recent progress, we still do not have 
precise knowledge of the constitutive law characterizing 
the stick-slip dynamics of earthquake faults. In fact, 
law of rock friction is often quite complicated, depending 
not just on the velocity or the displacement, but on the 
previous history and the "state" of contact surface, etc. 
The rate-and-state friction (RSF) law currently occupies 
the standard position among friction laws in the field 
of tectonophysics. Although the RSF law is formulated 
empirically three decades ago to account for certain as- 
pects of rock friction experiments (Dieterich, 1979; Ruina 
1983), the underlying physics was not known until very 
recently. While the RSF law shows qualitatively good 
agreement with numerous experiments, it is only good at 
aseismic slip velocities (slower than mm/sec). 

Among some progress made recently in the study of 
friction process, the most fascinating findings might be 
the rich variety of mechano-chemical phenomena, which 
comes into play at seismic slip velocities. Another impor- 
tant progress might be the understanding of the friction 
law of granular matter. This is also a very important 
point in understanding the friction law of faults as they 
consist of fine rock powder that are ground up by the 
fault motion of the past. The investigations on friction 
phenomenon at seismic slip velocities is now a frontier 
in tectonophysics. The RSF law no longer applies to 
this regime, where many mechano-chemical phenomena 
have been observed in experiments. The most illustrat- 
ing examples are melting due to frictional heat, thermal 
decomposition of calcite, silica-gel lubrication and so on. 
There have not been any friction laws that can describe 
such varied class of phenomena, which significantly af- 
fect the nature of sliding friction. In this review article, 
we wish to review the recent development concerning the 
basic physics of friction and fracture. 

Statistical physical study of earthquakes is usually 
based on models of various levels of simplification. There 
are several advantages in employing simplified models in 
the study of earthquakes. First, it is straightforward 
in the model study to control various material param- 
eters as input parameters. A systematic field study of 
the material-parameter dependence of real earthquakes 



3 



meets serious difficulties, because it is difficult to get pre- 
cise knowledge of, or even to control, various material pa- 
rameters characterizing real earthquake faults. Second, 
since an earthquake is a large-scale natural phenomenon, 
it is intrinsically not "reproducible" . Furthermore, large 
earthquakes are rare, say, once in hundreds of years for 
a given fault. If some observations arc to be made for 
a given large event, it is often extremely difficult to see 
how universal it is and to put reliable error bars to the 
obtained data. In the model, on the other hand, it is 
often quite possible to put reliable error bars to the data 
under well controlled conditions, say, by performing ex- 
tensive computer simulations. An obvious disadvantage 
of the model study is that the model is not the reality 
in itself, and one has to be careful in elucidating what 
aspect of reality is taken into account or discarded in the 
model under study. 

While numerous earthquake models of various levels 
of simplifications have been studied in the past, one 
may classify them roughly into two categories: The 
first one is of the type possessing an equation of mo- 
tion describing its dynamics where the constitutive rela- 
tion can be incorporated as a form of "force". The so- 
called spring-block or the Burridge-Knopoff (BK) model, 
which is a discretizcd model consisting of an assembly 
of blocks coupled via elastic springs, belong to this cat- 
egory (Burridge and Knopoff, 1967). Continuum models 
also belong to this category (Tse and Rice, 1986; Rice, 
1993). The second category encompasses further simpli- 
fied statistical physical models, coupled-lattice models, 
most of which were originally introduced as a model of 
SOC. This category includes the so-called Olami-Feder- 
Christensen (OFC) model (Olami, Feder and Chris- 
tensen, 1992), the fiber bundle model (Pradhan, Hansen 
and Chakrabarti, 2010), and the two fractal overlap 
model jBlrattacharva" et all 120091: iBhattacharvval 120051: 



IChakrabarti and Stinchcombd . 19991 ). These models pos- 
sesses extremely simplified evolution rule, instead of real- 
istic dynamics and constitutive relation. Yet, one expects 
that its simplicity enables one to perform exact or pre- 
cise analysis which might be useful in extracting essential 
qualitative features of the phenomena. 

It often happens in practice that, even when the 
adopted model looks simple in its appearance, it is 
still highly nontrivial to reveal its statistical properties. 
Then, the strategy in examining the model properties is 
often to perform numerical computer simulations on the 
model, together with the analytical treatment. In this re- 
view article, we wish to review the recent developments 
concerning the properties of these models mainly studied 
in statistical physics. 

Earthquake forecast is an ultimate goal of any earth- 
quake study. A crucially important ingredient playing 
a central role there might be various kinds of precur- 
sory phenomena. We wish to touch upon the follow- 
ing two types of precursory phenomena in this review 
article: The first type is a possible change in statisti- 
cal properties of earthquakes which might occur prior to 



mainshocks. The form of certain spatiotemporal corre- 
lations of earthquakes might change due to the proxim- 
ity effect of the upcoming mainshock. For example, it 
has been pointed out that the power-law exponent de- 
scribing the GR law might change before the mainshock, 
or a doughnut-likc quiescence phenomenon might occur 
around the hypoccntcr of the upcoming mainshock, etc. 
The second type of precursory phenomena is a possible 
nucleation process which might occur preceding main- 
shocks (Dieterich, 2009). Namely, prior to seismic rup- 
ture of a mainshock, the fault might exhibit a slow rup- 
ture process localized to a compact "seed" area, with its 
rupture velocity orders of magnitude slower than the seis- 
mic wave velocity. The fault spends a very long time in 
this nucleation process, and then at some stage, exhibits 
a rapid acceleration process accompanied by a rapid ex- 
pansion of the rupture zone, finally getting into a final 
seismic rupture of a mainshock. These possible precur- 
sory phenomena preceding mainshocks are of paramount 
importance in their own right as well as in possible con- 
nection to earthquake forecast. We note that similar nu- 
cleation process is ubiquitously observed in various types 
of failure processes in material science and in engineering. 

The purpose of the present review article is to help 
researchers link different branches of earthquake studies. 
First, we wish to link basic physics of friction and fracture 
underlying earthquake phenomenon to macroscopic prop- 
erties of earthquakes as a large-scale dynamical instabil- 
ity. These two features should be inter-related as "input 
versus output" or as "microscopic versus macroscopic" 
relation, but its true connection is highly nontrivial and 
still remains largely unexplored. To understand an ap- 
propriate constitutive law describing an earthquake in- 
stability and to make a link between such constitutive re- 
lations and the macroscopic properties of earthquakes is 
crucially important in our understanding of earthquakes. 
Second, we wish to promote an interaction between sta- 
tistical physicists and seismologists. We believe that the 
cooperation of scientists in these two areas would be very 
effective, and in some sense, indispensable in our proper 
understanding of earthquakes. 

Recently, there has been some progress made by sta- 
tistical physicists in characterizing the statistical aspects 
of the earthquake phenomena. These efforts arc of course 
based on established literature in seismology, physics of 
fracture and friction. Also, there has been considerable 
fusion and migration of the scientists and the established 
knowledge bases between physics and seismology. In this 
article, we intend to review the present state of our under- 
standing regarding the dynamics of earthquakes and the 
statistical physical modelling of such phenomena, start- 
ing with the same for fracture and friction. 

The article is organized as follows. In section II, we 
deal with the basic physics of fracture and friction. After 
reviewing the classic Griffith theory of fracture in II. A, 
we review a theory of fracture as a dynamical phase tran- 
sition in II. B. Rate and state dependent friction (RSF) 
law is reviewed in II. C, while the recent development be- 
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yond the RSF law is discussed in II. D. Section II. E is 
devoted to some microscopic statistical mechanical the- 
ories of friction. In section III, we deal with statistical 
properties of the model of our first type which includes 
the spring-block Burridge-KnopofF model (III. A) and the 
continuum model (III.B). In III. A, we examine statisti- 
cal properties of earthquakes including precursory phe- 
nomena with emphasis on both their critical and char- 
acteristic properties, while, in III.B, we mainly exam- 
ine characteristic properties of earthquakes and various 
slip behaviors including slow earthquakes. Implications 
of RSF laws to earthquake physics are also discussed in 
this section (III.B). In section IV, we deal with statistical 
properties of our second type of models which include the 
OFC model (IV. A), the fiber bundle model (IV.B) and 
the two fractal overlap model (IV. C). We also provide a a 
Glossary of some interdisciplinary terms as an Appendix. 



II. FRACTURE AND FRICTION 

A. GrIfFith energy balance and brittle fracture strength of 
solids 

In a solid, stress (a) and strain (S) bear a linear relation 
in the Hookean region (small stress). Non-linearity ap- 
pears for further increase of stress, which finally ends in 
fracture or failure of the solid. In brittle solids, failure 
occurs immediately after the linear region. Hence lin- 
ear elastic theory can be applied to study this essentially 
non-linear and irreversible phenomena. 

The failure process has strong dependence on, among 
other things, the diso rder properties of the material 
( Caldarelli et all Il999l). Often, stress gets concentrated 
arouii d the d isorder (IChakrabarti and Benguiguil Il997l : 
iLawnl . I1993I: iPetri et all Il994l ) where microcracks are 
formed. The stress values at the notches and corners 
of the microcracks can be several times higher than the 
applied stress. Therefore the scaling properties of disor- 
der plays an important role in breakdown properties of 
solids. Although the disorder properties tell us about the 
location of instabilities, it does not tell us about when a 
microcrack propagates. For that detailed energy balance 
study is needed. 

Griffith in 1920, equating the released elastic energy 
(in an elastic continuum) with the energy of the sur- 
face newly created (as the crack grows), arrived at a 
quantitative criterion for the equilibrium extension of the 
microcrack already present within the stressed material 
((Bergm an and Stroud . .1992i) . The following analysis is 
valid effectively for two-dimensional stressed solids with 
a single pre-existing crack, as for example the case of a 
large plate with a small thickness. Extension to three- 
dimensional solids is straightforward. 

Let us assume a thin linear crack of length 21 in an infi- 
nite elastic continuum subjected to uniform tensile stress 
a perpendicular to the length of the crack (see Fig. [1]). 
Stress parallel to the crack does not affect the stabil- 
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FIG. 1: A portion of a plate (of thickness w) under 
tensile stress a (Model I loading) containing a linear 
crack of length 21. For a further growth of the crack 
length by 2dl, the elastic energy released from the 
annular region must be sufficient to provide the surface 
energy ATwdl (extra elastic energy must be released for 
finite velocity of crack propagation). 



ity of the crack and has not, therefore, been considered. 
Because of the crack (which can not support any stress 
field, at least on its surfaces) , the strain energy density of 
the stress field {a'^ /2Y; where Y represents the elasticity 
modulus) is perturbed in a region around the crack, hav- 
ing dimension of the length of the crack. We assume here 
this perturbed or stress-released region to have a circular 
cross-section with the crack length as the diameter. The 
exact geometry of this perturbed region is not important 
here, and it determines only an (unimpo rtant) numeri- 
cal factor in the Griffith formula (see e.g. iLawnl (ll993l )V 
Assuming for the purpose of illustration that half of the 
stress energy of the annular or cylindrical volume, having 
the internal radius I and outer radius I -\- dl and length w 
(perpendicular to the plane of the stress; here the width 
w of the plate is very small compared to the other di- 
mensions), to be released as the crack propagates by a 
length d/, one requires this released strain energy to be 
sufficient for providing the surface energy of the four new 
surfaces produced. This suggests 

^{<T^ /2Y){2nwldl) > r(4u;dO. 

Here F represents the surface energy density of the solid, 
measured by the extra energy required to create unit sur- 
face area within the bulk of the solid. 

We have assumed here, on average, half of the strain 
energy of the cylindrical region having a circular cross- 
section with diameter 21 to be released. If this fraction 
is different or the cross-section is different, it will change 
only some of the numerical factors, in which we are not 
very much interested here. Also, we assume here linear 
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elasticity up to the breaking point, as in the case of brittle 
materials. The equality holds when energy dissipation, as 
in the case of plastic deformation or for the propagation 
dynamics of the crack, does not occur. One then gets 



A 

fr 



A 



(1) 



for the critical stress at and above which the crack of 
length 21 starts propagating and a macroscopic fracture 
occurs. Here A is called the critical stress-intensity factor 
or the fracture toughness. 

In a three-dimensional solid containing a single elliptic 
disk-shaped planar crack parallel to the applied tensile 
stress direction, a straightforward extension of the above 
analysis suggests that the maximum stress concentration 
would occur at the two tips (at the two ends of the ma- 
jor axis) of the ellipse. The Griffith stress for the brittle 
fracture of the solid would therefore be determined by the 
same formula ([T]), with the crack length 21 replaced by the 
length of the major axis of the elliptic planar crack. Gen- 
erally, for any dimension therefore, if a crack of length I 
already exists in an infinite elastic continuum, subject to 
uniform tensile stress a perpendicular to the length of 
the crack, then for the onset of brittle fracture , Griffith 
equates (the differentials of) the elastic energy Ei with 
the surface energy Es'. 



El 



a 

2Y 



r = Es 



Tl 



(2) 



where Y represents the elastic modulus appropriate for 
the strain, T the surface energy density and d the dimen- 
sion. Equality holds when no energy dissipation (due to 
plasticity or crack propagation) occurs and one gets 



(3) 



for the breakdown stress at (and above) which the exist- 
ing crack of length / starts propagating and a macroscopic 
fracture occurs. It may also be noted that the above for- 
mula is valid in all dimensions (d > 2). 

This quasistatic pic tu re ca n be extended 
( Pradhan and Chakrabartil l2003al) to fatigue be- 
havior of crack propagation for a < Of. At any stress a 
less than tr/, the cracks (of length /q) can still nucleate 
for a further extension at any finite temperature ksT 
with a probability exp[— E'/fcsT] and consequently 
the sample fails within a failure time r given by 



where 



exp[-£;(/o)/fcsT] 



E{1^) = E, +Ei^ Vll - —ll 



(4) 



(5) 



is the crack (of length /q) nuclcation energy. One can 
therefore express r as 



r-exp[A(l- ^)], 



(6) 



where (the dimensionless parameter) A ^ l^aj / {YksT) 
and CTj is given by Eq. This immediately suggests 
that the failure time r grows exponentially for a < aj 
and approaches infinity if the stress a is much smaller 
than af when the temperature ksT is small, whereas r 
becomes vanis hingly small as the stress a exceeds af, see, 
e.g.. lSornett^ (|2004 ) and lPoliti efd\ (|2002[ ). 



For disordered sohds, let us model the sohd by a perco- 
lating system. For the occupied bond/site concentration 
p > Pc, the percolation threshold, the typical pre-existing 
cracks in the solid will have the dimension (/) of correla- 
tion length £ ~ /S.p~'^ and th e elastic strength Y ~ Ap-^" 
( Stauffer and Aharonvl . ll99"2h . Assuming that the surface 
energy density F sc ales as ^'^^ , with the backbo ne (frac- 
tal) dimension ds I Stauffer and Aharonvl [T992h . equat- 
ing El and E^ as in one gets (|f ) ~ i'^'' ■ This 
gives 



with 



Tf^-[T, + {d-dBy] 



(7) 



for the 'average' fracture strength of a disordered solid (of 
fixed value) as one approaches the percolation thresh- 
old. Careful extensions of such scaling relations ([7]) 
and rigorous bounds for Tf has been obtained and com - 
pared extensively in Chakrabarti and Benguiguil ( 19971 ): 
iHerrmann and Rouxl (|l990l ): ISahimil (|2003l ). 



Extreme statistics of the fracture stress 

The fracture strength af oi a disordered solid does not 
have self-averaging statistics; most probable and the av- 
erage af may not match because of the extreme nature 
of the statistics. This is because, the 'weakest point' of 
a solid determines the strength of the entire solid, not 
the average weak points! As we have modelled here, the 
statistics of clusters of defects are governed by the ran- 
dom percolation processes. 

We have also discussed, how the linear responses, like 
the elastic moduli of such random networks, can be ob- 
tained from the averages over the statistics of such clus- 
ters. This was possible because of the self-averaging 
property of such linear responses. This is because the 
elasticity of a random network is determined by all 
the 'parallel' connected material portions or paths, con- 
tributing their share in the net elasticity of the sample. 

The fracture or breakdown property of a disordered 
solid, however, is determined by only the weakest (often 
the longest) defect cluster or crack in the entire solid. 
Except for some indirect effects, most of the weaker or 
smaller defects or cracks in the solid do not determine the 
breakdown strength of the sample. The fracture or break- 
down statistics of a solid sample is therefore determined 
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essentially by the extreme statistics of the most danger- 
ous or weakest (largest) defect cluster or crack within the 
sample volume. 

We discuss now more formally the origin of this ex- 
treme statistics. Let us consider a solid of linear size 
i, containing n cracks within its volume. We assume 
that each of these cracks have a failure probability 
fi{<j), i = 1,2, ... ,n to fail or break (independently) un- 
der an applied stress a on the solid, and that the per- 
turbed or stress-released regions of each of these cracks 
are separate and do not overlap. If we denote the cu- 
mulative failure probabilit y of the entire sample, un- 
der stress a, by F{a) then (IChakrabarti and Benguigui 



119971: iRav and Chakrabartilll985D 



1-F{a) =[](l-/,(a)) ~cxp 

i=l 

= exp [-L^~gia)] 



(8) 



where g{cr) denotes the density of cracks within the sam- 
ple volume L'^ (coming from the sum over the en- 
tire volume), which starts propagating at and above the 
stress level cr. The above equation comes from the fact 
that the sample survives if each of the cracks within the 
volume survives. This is the essential origin of the above 
extreme statistical nature of the failure probability F{a) 
of the sample. 

Noting that the pair correlation g{l) of two occu- 
pied sites at distance I on a percolation cluster de- 
cays as exp (— Z/^(p)), and connecting the stress cr with 
the length / by using Griffith's law (Eq. (P)) that 

a jk, one gets g{a) ^ exp(^-^^) for p p^. 

On sub stituting this, Eq. ([HI gives the G umbel distri- 
bution ( Chakrabarti and Beng uigui 1199/1) . If, on the 
other hand, one assumes a power law decay of g{l): 
g{l) ^ l^^ , then using the Griffith's law ([Ij, one gets 
g{a) ^ (^) , giving the Weibull distribution, from 
eqn. ([51), where m = b/a gives the Weibull modu- 
lus (jChakrabarti and Benguigu 1, [1993). The variation 
of F{a) with a in both the cases have the generic form 
shown in Fig.[2l F{a) is non-zero for any stress cr > and 
its value (at any a) is higher for larger volume {L''-). This 
is because, the possibility of a larger defect (due to fluc- 
tuation) is higher in a larger volume and consequently, 
its failure probability is higher. Assuming F[af) is finite 
for failure, the most probable failure stress af becomes a 
decreasing function of volume if extreme statistics is at 
work. 

The precise ranges of the validity of the Weibull or 
Gumbel distributions for the breakdown strength of dis- 
ordered solids are not well established yet. However, 
analysis of the results of detailed experimental and nu- 
merical studies of breakdown in disordered solids seem 
to suggest that the fluctuati ons of the extreme statistics 
dominate for small disorder ( Herrmann and Rouxl . [l990t 
ISahimil . [2OO3I) . Very near to the percolation point, the 
percolation statistics takes over and the statistics become 




FIG. 2: Schematic variation of failure probability F{(j) 
with stress cr for a disordered solid with volume Lf or 
Li {L2 > Li). 



self-a veraging. One can argue ( Bergman and Stroudl . 
Il992[ ). that arbitrarily close to the percolation thresh- 
old, the fluctuations of the extreme statistics will proba- 
bly get suppressed and the percolation statistics should 
take over and the most probable breaking stress becomes 
independent of the sample volume (its variation with dis- 
order being determined, as in Eqn. ([7]), by an appropriate 
breakdown exponent). This is because the appropriate 
competing length scales for the two kinds of statistics 
are the Lifshitz scale In L (coming from the finiteness of 
the volume integral of the defect probability: L'^{1 — pY 
finite, giving the typical defect size / ^ Ini) and the per- 
colation correlation length ^. When ^ < Ini, the above 
scenario of extreme statistics should be observed. For 
^ > InL, the percolation statistics is expected to domi- 
nate. 



B. Fracture as dynamical phase transition 

When a material is stressed, according to the linear elas- 
tic theory discussed above, it develops a proportional 
amount of strain. Beyond a threshold, cracks appear and 
on further application of stress, the material is fractured 
as it breaks into pieces. In a disordered solid, however, 
the advancing cracks may be stopped or pinned by the 
defect centers present within the material. So a compe- 
tition develops between the pinning force due to disorder 
and the external force. Upto a critical value of the ex- 
ternal force, the average velocity of the crack-front will 
disappear in the long time limit, i.e., the crack will be 
pinned. However, if the external force crosses this crit- 
ical value, the crack front moves with a finite velocity. 
This depinning transition can be viewed as a dynami- 
cal critical phenomena in the sense that near the crit- 
icallity universal scaling is observed which are indepen- 
dent of the microscopic d e tails o f the materials concerned 
(iBonamv and Bouchaudl . I2OIOI ). The order parameter 
for this transition is the average velocity v of the crack 
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Go " To 

FIG. 3: The average velocity of the crack front is 
plotted against external force {f = G — T, where G is 
the mechanical energy release rate and F is the fracture 
energy). For T ~ the depinning transition is seen. For 
finite temperature sub-crit i cal cr eep is shown 
( Ponson and Bonamv ^ 201Cl[). F rom 
( Ponson and Bonamvr2010l ). 



front. When the external force Z*^^* approaches the crit- 
ical value /^^* from a higher value, the order parameter 
vanishes as 



(r 



ext\6 



fr) 



(9) 



where 6 denotes the velocity exponent. It is to be men- 
tioned here that the pinning of a crack-front by disorder 
potential can occur at zero temperature. At finite tem- 
perature, there can be healing of cracks due to diffusion 
or there can be sub-c ritical crack propagation (in th e so 
called creep regime) (jBonamv and Bouchaudl [20101 ). In 
the later case, the velocity is expected to scale as 



exp 



-C 



ext \ 't> 



f 



(10) 



This sub -critica l scaling agrees well with experiments 
(jkoivisto et all [2007 : Ponson, 200^). In Fig. 11 the ex- 
perimental resul t of the crack propagation in the Botu- 
catu sandstone (jPonsonl . l2Q09| ) is shown. The average 
velocity of the crack is plotted against the mechanical 
energy release rate G (/ = G — F, where F is the fracture 
energy). The subcritical creep regime and the supercrit- 
ical power- law variations are clearly seen (insets) , which 
gives the velocity exponent close to sa 0.81 . 

Theoretical predictions of this exponent using Func- 
tional Renormalisations Group methods hav e placed its 
value around 9 = 0.59 ( Chauve et aLl . [200lh . where the 
experimental findings differ significantly {9 sa 0.80 ± 
0.15). Here we mention the numerical study of a model 
of the elastic crack-front propagation in a disordered 
solid. The basic idea is to consider the propagation of 
the crack front as an elastic string driven through a ran- 
dom medium. The crack front is characterised by an 




G (J.m- 



FIG. 4: Variation of average crack-front velocity against 
the mechanical en ergy release ra te is shown for 
Botucatu sandstone ( Ponsonl . |2009| ). The sub-critical 
creep region and supercritical power-law variations are 
shown in top-left and bottom right insets respectively. 
For the sub-critical regime, the data is fitted with a 
function v ~ g-'^/i^^-i'^))'^ for /x = 0.60 and 
(F) = Q5Jm~^. For the power-law variation 
(u ^ (G — GcY) in the super-cri tical region, 
Gc = UOJm~^ and 9 = 0.80. From (jPonsoni . [20091) . 



array of integral height (measured in the direction of the 
crack propagation) {hi, ft,2, . . . , /il} with periodic bound- 
ary conditions, where the unique values for the height 
profile suggest that any overhangs in the height profile is 
neglected. The forces acting on a site can be written as 



Mt) ^ ff + r' + gv^ih) 



(11) 



where f is the elastic force due to stretching, f^^* is 
the applied external force and 77 is due to disorder. The 
dynamics of the driven elastic chain is then given by the 
simple rule 



h,{t + I) ~ h^it) ^ v^it) 



1 if Mt)>Q 

otherwise. (12) 



The elastic force may have different forms in vari- 
ous contexts. When this force is short ranged (near- 
est neighbours) the well studied models are Edwards 



-Wilkinson (EW") (lEdwards and Wilkinson 


, Il98d) (see 


alsolAmar and Familv (1990l): Csahok et al. 


(119931)) and 


KPZ models (Kardar et all 11986^ Moser et al. Il991 


Sasamoto and SDohnl. 2010l) (see (Barabasi and Stanlev. 


1995) for extensive analysis). The lone 


range ver- 



sions includes the ones where the force decays as h > 
verse square (see e.g., (jPuemmer and Krauthl . l2007t )). 
The velocity exponent 9 , as is defined before , turn s 
out to be 0.625 ± 0.005 (jPuemmer and Krauthl . l2007t) . 
Also, mean field models (infinite range) are s t udied 
in this context (jLeschhornl . Il992l : IVannimenul [20021: 
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IVannimenus and Derridal l200ll ) (with 9 = 1/2; exactly). 
An infinite range model, where the elastic force only de- 
pends upon the total s tretching of the string, has als o 
been studied recently ( Biswas and Chakrabartil l201l[) . 
where the observed velocity exponent value {9 = 0.83 ± 
0.01) is r ather close to that found in some experiments 
(lPonsonl . [2009l) . 



C. Rate- and state-dependent friction law 

1. General remarks 

In a simplified view, an earthquake may be regarded as 
the rubbing of a fault. From this standpoint, friction 
laws of faults play a vital role in understanding and pre- 
dicting the earthquake dynamics. In addition, it should 
be noted that the celebrated Coulomb-Mohr criterion for 
brittle fracture involves the (internal) friction coefficient 
and thus the role of a friction law in earthquake physics is 
considerable. In this section, the phenomenology of fric- 
tion and its underlying physical processes are briefly re- 
viewed focusing on the recent developments. Some recent 
remarkable progress in experiments shall be also intro- 
duced, whereas, unfortunately, theoretical understand- 
ing of such experiments is rather poor. Thus, we try to 
propose the problems to be solved by physicists. 

Before explaining the knowledge obtained in the 20th 
and 21st centuries, it is instructive to see the ancient 
(16-17th centuries) phenomenology, which has been re- 
ferred to as the Coulomb- Amonton's law: (i) Frictional 
force is independent of the apparent area of contact, (ii) 
Frictional force is proportional to the normal load, (iii) 
Kinetic friction does not depend on the sliding velocity 
and is smaller than static friction. Among these three, 
the first two laws do not need any modification to this 
date, whereas the third law is to be modified and to be 
replaced by the celebrated rate- and state-dependent fric- 
tion law, which shall be introduced in the following sec- 
tions. The Coulomb- Amonton law in its original form 
is just a phenomenology involving only the macroscopic 
quantities such as apparent contact area and normal load. 
It is generally instructive to consider the sub-level (or mi- 
croscopic) ingredients that underlie such a macroscopic 
phenomenology. 

The essential microscopic ingredient in friction is as- 
perity, which is a junction of protrusions of t he surfaces 
(|Bowden and Tabod . l200ll: iRabinowiczl Il965l ). In other 
words, the two macroscopic surfaces in contact are in- 
deed detached almost everywhere except for asperities. 
The total area of asperities defines the true contact area, 
which is generally much smaller than the apparent con- 
tact area. Thus, the macroscopic frictional behavior is 
mainly determined by the rheological properties of as- 
perity. We write the area of asperity i as Ai. Then the 
total area of true contact reads 



A 



true — / ^ - 



(13) 



where S denotes the set of the asperities. This set de- 
pends on the state of the surfaces such as topography, 
and is essentially time dependent because the state of 
the surface is dynamic due to sliding and frictional heal- 
ing. 

Due to the stress concentration at asperities, molecules 
or atoms are directly pushed into contact so that an 
asperity may be viewed as a grai n boundary possibly 
with some inclusions and impurities ( Bowden and Taboii . 
l200ll: lRab inowiczl Il965() . Suppose that each asperity has 
its own shear strength ai, above which the asperity un- 
dergoes sliding. It may depend on the degree of grain- 
boundary misorientation and on the amount of impuri- 
ties at asperity. For simplicity, however, here we assume 

~ cry; i-c, the yield stress or shear strength of each 
asperity is the same. Then the frictional force needed to 
slide the surface reads 



F : 



ie5 



AiCTi 



(TY A, 



ay At 



(14) 



The frictional force is thus proportional to the area of 
true contact. Dividing Eq. ((T4)) by the normal force 
N, one obtains the friction coefficient /i = F/N. Using 
N = Ag,P, where A^, is the apparent area of contact and 
P is the normal pressure, one gets 



E 



Ai (Ji ^trucCy 



A^ P 



A,P 



(15) 



Alternatively, one can have Atruc/^a — fJ-P/'^Y- This 
means that the fraction of true contact area is propor- 
tional to the pressure normalized by the yield stress, 
where the friction coefficient is the proportionality co- 
efficient. Assuming that the yield stress of asperity is 
the same as that of the bulk, we may set ay ~ O.OIG, 
where G is the shear modulus. Inserting this and /i ~ 0.6 
into Eq. ([T5]), one has ArueMa ~ 60P/G. This rough 
estimate can be confirmed in experiment and numeri- 
cal si mulation ( Dieterich and Kilgorelll996t iHvun et all 
l2004l) . where the proportionality coefficient is on the or- 
der of 10. For example, at the normal pressure on the 
order of kPa, the fraction of true contact is as small as 
10-^ 



In view of Eq. (|T4)) . the first two laws of Coulomb- 
Amonton can be recast in the form that frictional force 
is proportional to the true contact area, which is inde- 
pendent of the apparent contact area but proportional to 
the normal load. This constitutes the starting point of 
a theory on friction, which shall be discussed in the fol- 
lowing subsections. The third law of Coulomb-Amonton 
is just a crude approximation of what we know of today. 
This should be replaced by the modern law, which is 
now referred to as the rate- and state- dependent friction 
(RSF) law. In the next subsection, we discuss the RSF 
law based on the first two laws of Coulomb-Amonton. 
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2. Formulation 

Extensive experiments on rock friction have been con- 
ducted in 1970s and 80s in the context of earthquake 
physics. An excellent review on these experimental works 
is done by Marone (1998). Importantly, these experi- 
ments reveal that kinetic friction is indeed not indepen- 
dent of sliding velocity. Thus, the third law of Coulomb- 
Amonton must be modified. Dieterich devises an empir- 
ical law that describes the behavior of friction coefficient 
(for both steady state and t ransient state) b ased on his 
experiments on rock friction ( Dietericbl . 1 1 979f ) . Later, the 
formulation is to some extent modified by Ruina (1983) 
by introducing additional variable(s) other than the slid- 
ing velocity. A new set of variables describes the state 
of the frictional surfaces so that they are referred to as 
the state variables. Although in general state variable(s) 
may be a set of scalars, in most cases a single variable is 
enough for the purpose. Hereafter the state variable is 
denoted by 9'. Using the state variable, the friction law 
reads 

f, = cJ + a'log^+b'\og^, (16) 

where a' and &' arc positive nondimcnsional constants, 
c' is a reference friction coefficient at a reference sliding 
velocity V^, , and £ is a characteristic length scale inter- 
preted to be comparable to a typical asperity length. In 
typical experiments, a' and b' are on the order of 0.01, 
and C is of the order of micrometers. Note that the state 
variable 9' has a dimension of time. 

The state variable 9' is in general time-dependent so 
that one must have a time evolution law for 9' together 
with Eq. (jlSp . Many empirical laws have been proposed 
so far in order to describe time-dependent properties of 
friction coefficie nt. One of th e commonly- used equations 
is the following (|Ruinal . [l983l) . 



(17) 



which is now referred to as the Dieterich's law or the 
aging law. This describes the time-dependent increase of 
the state variable even a.tV = 0. Meanwhile, other forms 
of evolution law may also be possible due to the empirical 
nature of Eq. p^ . For example, the follo wing one is als o 
known to be consistent with experiments ( Ruinal . [19831 ). 



V9' , V9' 



(18) 



which is referred to as the Ruina's law or the slip law. 
In a similar manner, a number of other evolution laws 
have been proposed so far, such as the composite of the 
slowness law and the slip law (Kato and Tulhs, 2001). 

Although there have been many attempts to clarify 
which evolution law is the most suitable, no decisive con- 
clusions have been made. As most of them give the iden- 
tical result if linearized around steady-sliding state, the 



difference between each evolution law becomes apparent 
only at far from steady-sliding state. One may immedi- 
ately notice that in Eq. (fTSl) the state variable is time- 
independent at y = so that it is not very quantitative 
in describing friction processes to which the healing is 
relevant. On the other hand, Eq. ^TE\i can describe a re- 
laxation process after the instantaneous velocity switch 
{V = Vi to V2) better than Eq. (H?]), while Eq. ^ 
predicts different responses for V — Vi to V2 and for 
= V2 to Vi , respectively. (Experimental data suggests 
that they are symmetric.) Also, it is known that Eq. 
(HH]) can describe a nuclea t ion p rocess better than Eq. 
pT]) (|Ampuero and Rubinl . l2008l ). However, we would 
not go further into the details of the experimental val- 
idation of evolution laws and leave it to the review by 
Marone (1998). 

Irrespective of the choice of evolution law, a steady 
state is characterized by 9'^^ = C/V so that the steady- 
state friction coefficient at sliding velocity V reads 



Mss = c' + (a' - 6')log-^ 



(19) 



Note that, as the nondimensional constants a' and b' are 
typically on the order of 0.01, the velocity dependence 
of steady-state friction is very small; change in sliding 
velocity by one order of magnitude results in ~ 0.01 (or 
even less) change in friction coefficient. It is thus natural 
that people in the 17th century overlooked this rather 
minor velocity dependence. However, this velocity de- 
pendence is indeed not minor at all but very important 
to the sliding instability problem: e.g., earthquakes. 

We also remark that Eq. p6|l together with an evo- 
lution law such as Eqs. (|17p or (fT8|). well describe the 
behavior of friction coe fficient not only for rock surfaces 
but als o metal surfa c es (IPodov et al\ . l201(jt) . two sheets of 
paper ( Heslot et al\ . \l99^ . etc. In this sense, the frame- 
work of Eq. ([T6|) is rather universal. This universality is 
partially because the deformation of asperities involves 
atomistic processes (i.e., creep). One can assume for 
creep of asperities ay = ksT /Ulog{V/Vo), where is 
an activation volume and Vq is a characteristic velocity 
involving the activation energy. Then Eq. leads to 



^ -ft rue 



\og{V/Vo) 



(20) 



where Ptruo = ^'^a/^tiuc is the actual pressure acting on 
the asperities. Comparing Eqs. ([20)) and with b' = 
(no healing) , one can infer that a' = J^p^ , as previously 



derived by some autho rs ( Heslot et al\ . 19941 : iNakatanl . 
I2OOII : iRice et all l200lD . However, we are unaware of a 
microscopic expression for b' to this date. We are also 
unaware of any microscopic derivations of evolution laws, 
such as Eqs. (|T7| . and (fTS)) . whereas an interesting effort 
to understan d the physical m eaning of evolution laws can 
be found in (lYoshiokaLfl997l) . 
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3. Stability of a steady state within the framework of the RSF 

As we discuss the earthquake dynamics based on the 
RSF, it is essential to discuss the frictional instabihty 
within the framework of the RSF. For simphcity, we con- 
sider a body on the frictional surface. The body is pulled 
by a spring at a constant velocity V. 



MX 



-k{X - Vt) - fiN, 



(21) 



where X is the position, M is the mass, k is the spring 
constant, and N is the normal load. This may be re- 
garded as the simplest model of frictional instability 
driven by tectonic loading. Suppose that the friction co- 
efficient fi is given by the RSF law Eq. together 
with an evolution law. The choice of the evolution law, 
i.e., Dieterich's or Ruina's, does not affect the follow- 
ing discussions as they are identical if linearized around 
steady state. The motion of the block is uniform in 
time if the surface is steady-state velocity strengthen- 
ing (a' > b') or if the spring constant is sufficiently large. 
For a steady-state velocity weakening surface the steady- 
sliding state undergoes Hopf bifurcation b elow a critica l 
spring constant. A linear stability analysis ( Heslot et all . 



spring constant. A Imear stability analysis l|Meslot et al.\ . 
ll994l : lR"uinal . Il983l ) shows that the steady-sliding state is 
unstable if 



(22) 



This relation plays a central role in various earthquake 
models, in which a constitutive law is given by the RSF 
law. This shall be discussed in section III. An important 
consequence of Eq. is that the tectonic motion is 

essentially stable ii a' — b' > 0. Namely, steady sliding 
is realized in the region where a' — b' > 0, whereas the 
motion may be unstable if a'^b' < 0. In addition, smaller 
£ widens the parameter range of unstable motion. 

Although the above analyses involve a one-body sys- 
tem, the stability condition Eq. appears to be es- 
sentially the same in many-body systems and continuum 
systems. Thus, provided that Eq. (PT|) applies, it is 
widely recognized in seismology that seismogenic zone 
has negative a b and smaller L, whereas aseismic zone 
has the opposite tendency. 



D. Beyond the RSF law 

It should be remarked that the RSF law has cer- 
tain limit of its application. Many experiments re- 
veal that the RSF no longer holds at high sliding 
velocities. This may be due to various mechano- 
chemical reactions that are induced by the frictional 
heat, which typically lubricate surfaces to a consid- 
erable degree; the friction coeffi cient becomes a s low 
0.2 or even less than 0.1 (iDi Toro et al 



as 



2004; 



Goldsbv and Tullisl. 120021: iHirose and Shimamotol. 12005 
Mizoguchi et al . 20061 : Tsutsumi and Shimamoto . ll99lt) . 
If such lubrication occurs in a fault, the fault motion is 



accelerated to a considerable degree and thus such effects 
have been paid much attention to during the last decade. 
Feedback of frictional heat may be indeed very important 
to faults, because the normal pressure in a seismogenic 
zone is of the order of 100 MPa. (Note that, however, the 
presence of high-pressure pore fluid may reduce the ef- 
fective pressure.) As this area of study is relatively new, 
the current status of our understanding on such mechano- 
chemical effects is rather incomplete. Taking the rapid 
development of this area into account, here we wish to 
mention some of the important experiments briefly. 



1. Flash heating 

Friction under such high pressure may lead to melting of 
rock. There have been some reports of molten rock ob- 
served in fault zones, which implies that the temperature 
is elevated up to 2000 K during earthquakes. 

A series of pioneering works on frictional melting in 
the context of earthquake has been conducted by Shi- 
mamoto and his coworkers. They devised a facility for 
rock friction at high speed under high pressure to find a 
behavior very different from that of the RSF law. Steady- 
state friction coefficient typically shows remarkable neg- 
ative dependence on slidi ng velocity and the relaxation 
to steady state is twofold ( Hirose and Shimamotol 120051 : 



iTsutsumi and Shimamo"tol . 19971) . At higher sliding ve- 
locity (e.g., 1 m/sec), friction coefficient decreases as low 
as 0.2 (or even less), whereas the typical value at qua- 
sistatic regime is around 0.7. We wish to stress that 
such a drastic decrease of friction coefficient cannot be 
explained in terms of the RSF law, where the change of 
steady-state friction coefficient is of the order of 0.01 even 
if the sliding velocity changes by a few orders of magni- 
tude, (recall Eq. p9)) . where a' and h' are both on the 
order of 0.01.) Thus, the mechanism of weakening must 
be qualitatively different from that of the RSF law. In- 
deed, in such experiments, molten rock is produced on 
surfaces due to the frictional heat. It is considered that 
so produced melt lubricates the surfaces to result in un- 
usually low friction coefficient. 

In view of Eq. ([T4)) , frictional melting must take place 
at asperities, where the frictional heat is produced. Thus, 
before the entire surface melts, asperities experience very 
high temperature, which may change the constitutive 
law. Such asperity heating has been also known in tribol- 
ogy and is referred to as flash heating. Rice applied this 
idea to fault friction in order to estimate the feasibility 
of flash he ating i n eart hquake dynamics. His argument is 
as follows (lRicel . l2006l) : The power input to asperity i is 
ayAiV, which is to be stored in the proximity of asperity. 
As discussed later, it is essential to assume here that heat 
conduction is one-dimensional; i.e., temperature gradi- 
ent is normal to the surface, whereas uniform along the 
transverse directions. The produced heat invades toward 
the bulk over the distance yi^tht, where i?th is thermal 
diffusivity. Thus, frictional heat is stored in the small 
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volume of Ai^Jat. Writing the average temperature of 
this hot volume as T'(i), the deposited thermal energy 
reads cppAi^/Dtht{T{t) — Tq), where cp is the isobaric 
specific heat, p is the mass density, and Tq is the ambient 
temperature. Then the energy balance leads to 



T{t) - To 



(TyV / t 



pep 



D 



th 



(23) 



This indicates that the surface temperature increases 
with time as ^/t. Writing as the critical tempera- 
ture above which an asperity looses its shear strength, 
then the duration for the temperature to be elevated 
up to the critical temperature reads 



th 



PCpjTyj - Tq) 

ctyV 



(24) 



This heating process is limited to the duration or life- 
time of an asperity contact. If we write the longitudinal 
dimension of each asperity as Ci , the lifetime of an asper- 
ity is estimated as Ci/V . Thus, weakening of an asperity 
occurs if and only if < Ci/V. Taking Eq. (p4|) into 
account, this condition may be written as 



V > 



D 



th 



PCp{Tyj - Tq) 



<7Y 



(25) 



Neglecting the statistics of £i, one gets the characteristic 
sliding velocity 14, above which weakening occurs. 
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C 



pCp{Tw - Tq) 



ay 



(26) 



Alternatively, from Eq. (j25p . the maximum size of 
asperity that does not melt at the sliding velocity V is 
given by 



Ah 
V 



pcp{Tw - Tq) 



cry 



(27) 



The proportion of non-melting asperity may be approxi- 
mated by Lniax/C Assuming that the friction coefficients 
of a molten asperity and a non-melting one are given as 



r/i, {T<T^ 

\/2. (r>T,„: 

the average friction coefficient reads 



C 



/2 (^1 



/2 + (/l - /2) 



V 



(28) 

(29) 
(30) 



The friction coefficient decreases a.s V ^ a.t high slip ve- 
locity V >Vio. Taking a = Imm^/s, pcp = 4 MJ/m^K, 
D = 5/im, r„ - To = TOOK, and ay = 0.02 to 0.1 G 
(shear modulus) = 0.6 to 3 GPa, the characteristic ve- 
locity Vw is 0.5 to 14 m/s. This does not contradict rock 



experiments on melting-induced weakening. Also, com- 
parison of Eq. (j30p with experiments is not inconsistent, 
although /i and /2 are fitting parameters. 

Note that the above discussion does not depend on the 
apparent normal pressure, as the pressure on asperity is 
approximately the yield stress (of uniaxial compression) 
irrespective of the apparent normal pressure. Thus, fiash 
melting could occur in principle even when the appar- 
ent pressure is very low as long as the sliding velocity is 
larger than given by Eq. (j26p . However, in an experi- 
ment conducted at relatively low pressures, the threshold 
velocity is an order of magnitude smaller t han t he pre- 
diction of Eq. ^ (jKuwano and Hatanol l201l[) . This 
may be because other relevant mechanisms are responsi- 
ble for dynamic weakening observed in experiments, but 
the answer is yet to be given. 

It is also important to notice that in the above discus- 
sion the assumption of one dimensional heat conduction 
is essential; i.e., the frictional heat is not transferred in 
the horizontal directions but only to the normal direc- 
tion. This assumption implies that the thermal diffusion 
length \J cA^, must be smaller than the height of a pro- 
trusion that constitutes an asperity. Assuming that the 
height of a protrusion is proportional to a horizontal di- 
mension L, , this condition leads to 



DthPCpjT^, - Tq) 

(JyV 



(31) 



Because it is estimated in general that pcp{T.uj — Tq) > 
ay, Eq. ([5T|) immediately follows from Eq. (^5]) . Thus, 
the assumption of one-dimensional heat conduction may 
be sound. 

If the asperities are sufficiently small so that the ther- 
mal diffusion length exceeds the height of protrusions, 
the assumption of one-dimensional heat conduction is vi- 
olated. A good example is friction of nanopowders, in 
which a typical size o f the true c ontac t area is on the 
order of nanometers ( Han et oil . l201lf ). Interestingly, 
one can still observe dynamic weakening similar to those 
caused by ffash melting, but they did not attribute this 
behavior to ffash melting, because the duration of contact 
between nano-grains was too short to cause the signiff- 
cant temperature increase. The physical mechanism of 
such weakening is still not clear. (Silica-gel lubrication 
may be ruled out as the material they used is silica- free.) 



2. Frictional melting and thermal pressurization 

There is yet another class of weakening phenomena called 
frictional melting; the melt is squeezed out of asperities 
to ffU the aperture between the two surfaces. Such a situ- 
ation can occur if the surfaces are rubbed for sufficiently 
long time. If this process occurs, the melt layer supports 
the apparent normal pressure to reduce the effective pres- 
sure at asperities and ultimately hinders solid-solid con- 
tact. This leads to the disappearance of the asperities; 
i.e., no solid-solid contacts between the surfaces but a 
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thin layer of melt under shear. There are some analyses 
of such systems assuming the Arrhenius-type viscosity 
(jFialko and Khazanl . 12003: iNielsen et al ]. 120081) . [n do- 
ing so, one can predict the shear traction is proportional 
to P^/^, where P is the normal pressure. The quantita- 
tive validation of such theories is yet to be done. 

It may be noteworthy here that the viscosity of such a 
liquid film involves rather different problem; nanofluidics. 
The melt may be regarded as nanofluid, the viscosity of 
which may be very different from the ordinary ones. The 
shear flow of very thin layers of melt (under very high 
press ure) may be unstable d ue to the partial crystalliza- 
tion (|Thompson et all Il992h Till this date, the effect of 
nanofluidics on frictional melting is not taken into ac- 
count and left open to physicists. 

Meanwhile, the evidence of frictional melting of a fault 
is not very often found in core samples or in outcrops. As 
faults generally contain fluid, frictional heat increases the 
fluid temperature as well. As a result, the fluid pressure 
increases and the effective pressure on solid-solid contact 
decreases. Therefore, the frictional heat production gen- 
erally decreases in the presence of fluid. In the simplest 
cases where the fault zone is impermeable, the effective 
friction (and the produced frictional heat) may vanish as 
the fluid pres sure can be as large as the rock pressure 
(|Sibsonl . \l97?t\ . This is referred to as thermal pressur- 
ization, and a large number of work has been devoted to 
such dynamic interaction between frictional heat and the 
fluid pressure. More detailed formulations incorporate 
the effect of fluid diffusion wit h nonzero permeability o f 
host rocks (|Lachenbruchl Il980l: iMase and Smithl . Il987t) . 
In these analyses, the extent of weakening is enhanced if 
a fault zone has smaller compressibility and permeabil- 
ity. Although this behavior is rather trivial in a qual- 
itative viewpoint, some nontrivial behaviors are found 
in a model where the permeability is assumed to be a 
dynamic quantity coup l ed wi th the total displacement 
( Suzuki and Yamashital . r2010| ). However, it is generally 
difficult to judge the validity of a model from observa- 
tions and thus we do not discuss this problem further. 



3. Other mechanochemical effects 

In some systems, anomalous weakening of friction (p ^ 
0.2) can be observed at sliding velocities much lower than 
the critical velocity for flash heating (Eq. ([26])). Typi- 
cally, one can observe weakening at sliding velocity of 
the order of mm/s. Thus, there might be mechanisms 
for drastic weakening other than frictional melting. 

Such experiments are typically conducted with com- 
plex materials like fault gouge taken from a natural fault 
so that there may be many different mechanisms of weak- 
ening depending on the specific compositions of rock 
species. Among them, the mechanism that might bear 
some robustness is the lu brication by silica-gel pro duc- 
tion (|Di Toro et all |2004 iGoldsbv and Tulli i. l2002[ ). In 
several experiments on silica-rich rock such as granite. 



SEM observation of the surfaces reveals a silica-gel layer 
that experienced shear flow. The generation of silica-gel 
may be attributed to chemical reactions between silica 
and water in the environment. This silica-gel is inter- 
vened between the surfaces to result in the lubrication 
of fault. Although the details of the chemical reactions 
is not very clear, the mechanics of weakening may be 
essentially the same as that of flash heating and melt- 
ing, because in the both cases the cause of weakening 
is some soft materials (or liquids) that are produced by 
shear and intervened at asperities. However, in the case 
of silica-gel formation, the thixotropic nature of silica-gel 
may result in peculiar behaviors of friction, as observed 
in experiment by Di Toro et al (2004). 

In addition, we wish to add several other mechanisms 
that lead to anomalous weakening. Han et al. (2007) 
found friction coefficient as low as 0.06 in marble under 
relatively high pressure (1.1 to 13.4 MPa) and high slid- 
ing velocity (1.3 m/s). Despite the utilization of several 
techniques for microstructural observation, they could 
not observe any evidence for melting such as glass or 
amorphous texture but only a layer of nanoparticles pro- 
duced by thermal decomposition of calcite due to fric- 
tional heating. Mizoguchi et al. (2006) also found fric- 
tion coefficient as low as 0.2 in fault gouge taken from 
a natural fault, where they also could not find any evi- 
dences for melting. To this date, the mechanism of such 
frictional weakening at higher sliding velocity is not clear. 

It might be important to notice that these samples 
inevitably include a large amount of sub-mi cron grains 



that are worn out by high- s peed friction ( Han et oil . 
l2007l : iHavashi and Tsutsumil . l2010t) . They may play 
an important role in weakening at high sliding veloci- 
ties. The grain size distribution of fault gouge is typ- 
ically wel l fitted by a power l aw wi th exponent —2.6 
to —3.0 ( Chester and Chested . Il998t ) so that smaller 
grains cannot be neglected in terms of volume frac- 
tion. The exponent appear s to be common to laboratory 
( Marone and Schol4 Il989l) or numerical experiments of 



wear ( Abe and Mairl . |2009[ ). Rheology of such fractal 
grains has not been investigated in a systematic man- 
ner, notwithsta nding a pioneering computational work 
(lMorganl[l999l) . The influence of grains to friction shall 
be discussed in detail in the next subsection. 



4. Effect of the third body: granular friction 

Previously, we considered the situation where two sur- 
faces were in contact only at asperities. This is generally 
not the case if the asperities are worn out to be free parti- 
cles that are intervened between the two surfaces. In this 
case, a system can be regarded as granular matter that is 
sheared by the two surfaces. The core of a natural fault 
alway s consists of powdered rock ( Chester and Chesteil . 
119981) . which is produced by the fault motion of the past. 
Thus, friction on fault is closely related to the rheology 
of granular rock. 
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As is briefly mentioned in the previous subsection, 
earthquake physics involves a wide range of sUding ve- 
locities (or shear rate) ranging from tectonic time scale 
(e.g., nm/s) to coseismic scale (m/s). It is thus plau- 
sible that the rheological properties of granular matter 
is qualitatively different depending on the range of slid- 
ing velocities. Here we define two regimes for granu- 
lar friction: quasititatic and dynamic regimes. In the 
quasistatic regime, the frictional properties of granular 
matter is described by the RSF. However, some impor- 
tant properties will be remarked that are not observed for 
bare surfaces. In the dynamic regime, one may expect dy- 
namic strength e ning as observed i n num erical simulations 
( da Cruz et 

M 120051 : iGDRMidil . 120041 ). However, at the 
same time one may also expect weakening due to vari- 
ous mechanochemical reacti ons ( Havashi and Tsutsumil . 
I2OIOI : iMizoguchi et aLl . 120061 ). The rheological properties 
that are experimentally observed are determined by the 
competition of these two ingredients. Here we review 
the essential rheological properties of granular matter in 
these two regimes. 

In experiments on quasistatic deformation, friction of 
granular matter seems to obey the RSF law. However, 
some important properties that are different from those 
of bare surfaces should be remarked. 

1. Velocity dependence of steady-state friction ap- 
pears to be affected by the layer thickness. In 
particular, the value of a' — b' in Eq. (|19p is an 
increasing function of the layer thickness. 

2. The value of a' — b' appears to have negative de- 
pendence on the total displacement applied to a 
system. This is true both for granular matter and 
bare surfaces. 

3. Transient behaviors can be described by cither Di- 
etcrich's or Ruina's laws, as in the case of bare sur- 
faces. The characteristic length in an evolution law 
is proportional to the layer thickness. 

These experimental observations arc well summarized 
and discussed in detail by Marone (1998). We thus shall 
not repeat them here and just remark the essential points 
described above. 

As to the first point, there is no plausible explanation 
to this date. It appears that the second point could be 
merged into the first point if the effective layer thickness 
(i.e., the width of shear band) decreases as the displace- 
ment increases. However, we wish to remark that it is 
also true in the case of bare surfaces, where the effec- 
tive layer thickness is not a simple decreasing function 
of the displacement. Thus, the second point cannot be 
explained in terms of the thickness. The third point indi- 
cates that the shear strain is a more appropriate variable 
than the displacement of the boundary for the descrip- 
tion of the time evolution of friction coefficient. This 
may be reasonable as the duration of contact between 
grains is inversely proportional to the shear rate. How- 
ever, the derivation of evolution laws (either Dietcrich's 



or Ruina's) from the grain dynamics is not known to this 
date. To construct a theory that can explain these three 
properties based on the nature of granular matter is still 
a challenge to statistical physicists. 

Then we discuss the dynamic regime. Rheology of 
granular matter in the dynamic regime is exte n sive^ 



granular matter m tnc dynamic regime is exte n sively 
investigated in statistical physics ( da Cruz et all 120051 : 
[CDRMidi. 2004i) . As to the steady-state friction coefla- 
cient, the shear-rate dependence is one of the main in- 
terests also in statistical physics. There are many in- 
gredients that potentially affect the friction coefficient 
of granular matter: the grain shape, degree of inelastic- 
ity (coefficient of restitution), friction coefficient between 
grains, the stiffness, pore-fluid, etc. Shape dependence is 
very important to granular friction, but theoretical un- 
derstanding of this effect is still very poor. Thus, for sim- 
plicity, we neglect the shape effect and involve only spher- 
ical grains. Furthermore, we limit ourselves to the effects 
of shear rate, stiffness, mass and diameter of grains, coef- 
ficient of restitution, and intergrain friction. This means 
that we neglect time- or slip- depen dent deformation of 
the g rain contacts, such a s wear (iMarone and Schold . 
119891 ) or frictional healing (jBocauet et all Il998l) . The 
effects of pore-fluid is also neglected; i.e., we discuss only 
the dry granular matter here. With such idealization, one 
can make a general statement on a constitutive law by 
dimensional analysis. The friction coefficient of granular 
matter is formally written as 



fi = ^(P,m,d,7, r, /Xe,e), 



(32) 



where P is the normal pressure, m is the mass, d is the 
diameter, 7 is shear rate, Y is the Young's modulus of 
grains, fie is the intergrain friction coefficient, and e is 
the coefficient of restitution. (It should be noted that one 
assumes a single characteristic diameter d in Eq. ((32)) .) 
From the viewpoint of dimensional analysis, the argu- 
ments on the left hand side of Eq. ([5^ must be nondi- 
mensional numbers. 



M = /^(^,K,Me,e), 



(33) 



where / ~ jy^m/Pd and k = Y/P. Thus, the friction co- 
efficient of granular matter depends in principle on these 
four nondimensional parameters. Many numerical simu- 
lations reveal that fi is rather insensitive to k and e, and 
the shear-rate dependence is mainly described by /. This 
nondimensional number / is referred to as the inertial 
number. Importantly, t he dependence on / i s positive in 
numerical simulations ( da Cruz et Ml, 120051: iGDRMidi 
I2OO4I : iHatand . l2007t ) ; namely, the shear-rate dependence 
is positive. It is important to remark that the negative 
shear-rate dependence, which is ubiquitously observed in 
experiments, cannot be reproduced in numerical simu- 
lation. This is rather reasonable because the origin of 
the negative velocity dependence is the time dependent 
increase of true contact, whereas in simulation the pa- 
rameters are time- independent. 

Experiments in the context of earthquake physics are 
conducted at relatively high pressures at which the fric- 
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tional heat affects the ph ysical state of granular mat- 
ter. In some experime nts ( Havashi and Tsutsumil . I2OIOI : 
iMizoguchi et al\ . l2006l ) . remarkable weakening (fi ~ 0.1) 
is observed. Because such anomalous behaviors may in- 
volve shear-banding as well as various chemical reactions 
such as thermal decomposition or silica-gel formation, 
the frictional properties should depend on the detailed 
composition of rock species contained in granular matter. 
These weakening behaviors must be further investigated 
by extensive experiments. 

So far we discussed the steady-state friction coefficient, 
but the description of transient states are also important 
in understanding the frictional instability (and earth- 
quake dynamics) . An evolution law for quasistatic regime 
is indeed essentially the same as that for ba re surfaces; 
namely, aging law or slip law ( Maronel Il998l ) . An evolu- 
tion law in the dynamic regime is well described by the 
linear relaxation equation even for relatively large veloc- 
ity change (Ha tano. 200^. 



(34) 



where fiss is the steady-state friction coefficient that de- 
pends on the sliding velocity and r is the relaxation time. 
It is found that t is the relaxation time of the velocity 
profile inside granular matter and is scaled with ^JmJPd. 
Thus, importantly, the inertial number, which describes 
steady-state friction may be written using r as / ~ T7; 
i.e., the shear rate multiplied by the velocity relaxation 
time. It may be noteworthy that the inertial number is 
an example of Deborah number, which is in general the 
internal relaxation time normalized by the experimental 
time scale. 

Note the difference from the conventional evolution law 
in the framework of the RSF law; Eqs. ^\ and ([TH]). It 
is essential that Eq. (p4)) does not contain any length 
scale but only the time scale. This means that the relax- 
ation process of high-speed granular friction takes time 
rather than the slip distance. However, we wish to stress 
that the validity of Eq. (j34[) is found only in simulation 
on dry granular matter and is not verified in physical 
experiment. 



E. Microscopic theories of friction 

Many attempts were made in explaining friction from 
an atomistic point of view. Of course, such effort are 
meaningful only when the surfaces are smooth and the 
atomistic properties determine friction. This approach 
has gained importance in recent years, because of ad- 
vancement of technology in this field. Due to Atomic 
Force Microscopy (AFM) etc., sliding surfaces can now 
be probed upto atomic scales. Also, present day com- 
puters allow large scale molecular dynamics simulation 
that helps in understanding the atomic origin of fric- 
tion. In this approach, the ato mic origin of f r iction 
forces are investigated (see also (iBhushan et ai . Il995t 
iBraun and Naumovetsl . l2006l : iHolscher et all \20m ). In 



this purpose, two atomically smooth surfaces arc taken 
and by writing down the equations of motion, friction 
forces are calculated. Effect of inhomogeneity, impurity, 
lubrication and disorder in terms of vacancies of atoms 
arc also considered. 

One of the foremost attempts t o model friction from 
atomic origin was of Tomlinson's (iTomlinsonl . Il929l) . In 
this model, only one atomic layer of the surfaces in con- 
tact are considered. In particular, the lower surface in 
considered to be rigid and provide a periodic (sinusoidal) 
potential for the upper body. The contact layer of the 
upper body is modelled by mutually disconnected beads 
(atoms) which are attached elastically to the bulk above. 
This model is, of course, oversimplified. The main draw- 
back is that no interaction between the atoms of the up- 
per body is considered. 



1. Frenkel Kontorova model 



Frenkcl-Kontorova ( Frenkel and Kontoroval [19381) model 
overcomes some of these difficulties. In this model the 
surface of the sliding object is modelled by a chain of 
beads (atoms) connected harmonically by springs. The 
base is again represented by a sinusoidal potential. The 
Hamiltonian of the system can, therefore, be written as 



N 1 



if + Vix,)], (35) 



where, Xi is the position of the i-th. atom, a is the equi- 
librium spacing of the chain and V{x) = —Vocos{^^). 
Clearly, there arc two competing lengths in this model, 
viz. the equilibrium spacing of the upper chain (a) and 
the period of the substrate potential (6). While the first 
term tries to keep the atoms in their original positions, 
the second term tries to bring them in the local minima 
of the substrate potential. Simultaneous satisfaction of 
these two forces is possible when the ratio a/6 is com- 
mensurate. The chain is then always pinned to the sub- 
strate in the sense that a finite force is always required 
to initiate sliding. Below that force, average velocity 
vanishes at large time. However, interesting phenom- 
ena takes place when the ratio a/6 is incommensurate. 
In that case, upto a finite value of the amplitude of the 
substrate potential, the chain remains "free" . In that 
condition, for arbitrarily sm all external force, sliding is 
initiated. The hull function ( Pevrard and Aubrvl . 1 19831 ) 
remains analytic. Beyond the critical value of the ampli- 
tude, the hull function is no longer analytic and a finite 
external force is now required to initiate sliding. This 
transition is called the b reaking of analyticity trans ition 



transition is called tne b reaking or analyticity trans ition 
or the Aubry transi tion (iPevrard and Aubrvl. 1983| ) (for 
extensive details see iBraun and Kivshail (|2004[ )V 
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2. Two-chain model 

The Frenkel-Kontorova model has been generaUsed in 
many ways viz., extension in higher dimensions, effect 
of impurit y, the Frenkel-Kontorova - Toml inson model and 
so on (seeH raun and NaumovetsI (|2006f ) and references 
therein). But one major shortcoming of the Frenkel- 
Kontorova model is that the substrate or the surface 
atoms of the lower substance are considered to be rigidly 
fixed in their equilibrium position. But for the same rea- 
son why the upper surface atoms should relax, the lower 
surface atoms should relax too. In the two chain model of 
friction ( Matsukawa and Fukuvanial . Il994l ) this question 
is addressed. In this model, a harmonically connected 
chain of atoms is being pulled over another. The atoms 
have only one degree of freedom in the direction parallel 
to the external force. The equations of motion of the two 
chains are 



Ka{Xi+i + - 2Xi) 

+ ^Fi{x,~yj)+Fe,, (36) 



+ ^ Filjji - Xj) - Ksiui - icjST) 

where, Xi and yi denotes the equilibrium positions of 
the upper and lower chain respectively, to's represent 
the atomic masses, 7's represent the dissipation constant, 
K''s the strength of inter-atomic force and iV's the num- 
ber of atoms in each chain, c's the lattice spacing, while 
sufhx a denotes upper chain and suffix b denotes the lower 
chain. F^x is the external force and Fj is the inter-chain 
force between the atoms, which is derived from the fol- 
lowing potential 



2 Cb 



(38) 



where Kj is the interaction strength. 

It is argued that the frictional force is of the form 



(39) 



It is then shown by numerical analysis that the velocity 
dependence of the kinetic frictional force becomes weaker 
as the static friction increases (tuned by different K^s). 
The velocity dependence essentially vanishes when static 
frictional force is increased, giving one of the Amonton- 
Coulomb laws. 

In this case, the lower chain atoms, which forms the 
substrate potential, is no longer rigidly placed. Still, the 
breaking of analyticity transition is observed. Fig. [5] 
shows the variation of the maximum static frictional force 
with interaction potential strength. For different values 
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FIG. 5: The variations of the maximum static friction 
with the amplitude (Ki) of the inter-chain potential for 

different values of the lower-chain stiffness {Kg). The 
limit Kg — > 00 corresponds to Frenkel-Kontorova model. 
But it is clearly seen that even for finite Ks (i.e., when 
the lower chain can relax) there is a finite value of 

inter-chain potential amplitude upto which the static 

friction is practically zero and it increases afterwards, 
sig nifying Aubry transition (iMats ukawa and Fukuyamal 
119941) . From (jMatsukawa and Fukuvamal . ll9Ml 



of the rigidity with which the lower chain is bound (Ks), 
different curves are obtained. This indicates a pinned 
state even for finite rigidity of the lower chain. 



3. Effect of fractal disorder 

Effects of disorder and impurity have been studied in the 
microscopic models of friction. Also there have been ef- 
forts to incorpor ate the effect of self -affine roughness in 
friction. In Ref. ( Eriksen et al I I2OIOD . the effect of disor- 
der on static friction is considered. A two-chain version 
of the Tomlinson model is considered. The self-affine 
roughness is introduced by removing atoms and keep- 
ing the remaining ones arranged in the form of a Cantor 
set. The Cantor set, as is discussed before, is a simple 
prototype of fractals. Instead of considering the regular 
Cantor set, here a random version of it is used. A line 
segment [0,1] is taken. In each generation, it is divided 
into s equal segments and s — r of those are randomly re- 
moved. In this way, a self similar disorder is introduced, 
which is present only in the statistical sense, rather than 
strict geometric arrangement. 

This kind of roughness is introduced in both the chains. 
Then the inter-chain interaction is taken to be very short 
range type. Only when there is one atom exactly over the 
other (see Fig. [5]) there is an attractive interaction. In 
this way, the maximum static friction force can be cal- 
culated by estimating the overlap of these two chains. It 
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III. EARTHQUAKE MODELS AND STATISTICS I: 
BURRIDGE-KNOPOFF AND CONTINUUM MODELS 

In the previous section, we reviewed the basic physics of 
friction and fracture, which constitutes a "microscopic" 
basis for our study of "macroscopic" properties of an 
earthquake as a large-scale frictional instability. Some 
emphasis was put on the RSF law now regarded as the 
standard constitutive law in seismology. In this and fol- 
lowing sections, we wish to review the present status of 
our research on various types of statistical physical mod- 
els of earthquakes introduced to represent their "macro- 
scopic" properties. 



V(x) 



FIG. 6: Schematic representation of the two chain 
version of the Tomlinson model with (a) no disorder, 
(b) Cantor set diso rder, (c)the ef f ective substrate 
potential ( Eriksen et al\ . \201(i }. 




FIG. 7: The overlap distribution for s = 9, r = 8 is 
shown. The dotted curve shows the average distribution 

with random off-set and the continuous curve is that 
without random off set. The distribution is qualitatively 
different from the Gauss ian distribution expe cted for 
random disorder ( Eriksen et all |2010[ ) . 



turns out that the static friction force has a distribution, 
which is qualitatively different from what is expected if a 
random disorder or no disorder is present. The scaled (in- 
dependent of gene ration) distribution o f overlap or static 
friction looks like ( Eriksen et aLl . [2010[ ) 



r-^^{x/R)/R = J2 c'-'^imr''' l terms . . . rHx), 

(40) 



where R = j s and (f''^{x) 



For a par- 



ticular (s,r) combination (9,8), the distribution function 
is shown in Fig. [T] It clearly shows that the distribu- 
tion function is qualitatively different from the Gaussian 
distribution expected if the disorder were random. 



A. Statistical properties of the Burridge-Knopoff model 

1. The model 

One of the standard models widely employed in statisti- 
cal physical study of earthquakes might be the Burridge- 
Knopoff (BK) model (Rundle, 2003; Ben-Zion, 2008). 
The model was first introduced in Burridge and Knopoff 
(1967). Then, Carlson, Langer and collaborators per- 
formed a pioneering study of the statistical properties 
of the model (Carlson and Langer, 1989a; Carlson and 
Langer, 1989b; Carlson et al., 1991; Carlson, 1991a; 
Carlson, 1991b; Carlson, Langer and Shaw, 1994), pay- 
ing particular attention to the magnitude distribution of 
earthquake events and its dependence on the friction pa- 
rameter. 

In the BK model, an earthquake fault is simulated by 
an assembly of blocks, each of which is connected via the 
elastic springs to the neighboring blocks and to the mov- 
ing plate. Of course, the space discretization in the form 
of blocks is an approximation to the continuum crust, 
which could in principle give rise to an artificial effect 
not realized in the continuum. Indeed, such a criticism 
against the BK model employing a certain type of fric- 
tion law, e.^., the purely velocity- weakening friction law 
to be defined below in 2[A], was made in the past (Rice, 
1993), which we shall return to later. 

plate drive 

//////////^ TTT^ 




friction 



FIG. 8: The Burridge-Knopoff (BK) model in one 
dimension. 
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In the BK model, all blocks are assumed to be sub- 
jected to friction force, the source of nonlinearity in the 
model, which eventually realizes an earthquake-like fric- 
tional instability. As mentioned in section II, the stan- 
dard friction law in modern seismology might be the RSF 
law. In order to facilitate its computational efficiency, 
even simpler friction law has also been used in simula- 
tion studies made in the past. 

We first introduce the BK model in one dimension 
(ID). Extension to two dimensions (2D) is straightfor- 
ward. The ID BK model consists of a ID array of N 
identical blocks, which arc mutually connected with the 
two neighboring blocks via the clastic springs of the elas- 
tic constant kc, and are also connected to the moving 
plate via the springs of the elastic constant fcp, and are 
driven with a constant rate: See Fig. [8] All blocks are 
subjected to the friction force $, which is the only source 
of nonlinearity in the model. The equation of motion for 
the i-th block can be written as 

mC/, ^ kp{v't' - U,) + fc,([/,+i - 2U, + U,-i) - (41) 

where f is the time, Ui is the displacement of the i-th 
block, I'' is the loading rate representing the speed of 
the moving plate, and $i is the friction force at the i-th 
block. 

In order to make the equation dimensionless, the time 
t' is measured in units of the characteristic frequency uj = 
y^kp/m and the displacement Ui in units of the length 
L* ~ ^o/kp, $0 being a reference value of the friction 
force. Then, the equation of motion can be written in 
the dimensionless form as 



ut ~ Ui + P{ui-f-i — 2ui 4- Ui^i) 



(42) 



where t = t'ut is the dimensionless time, Ui = Ui/ L* is the 
dimensionless displacement of the i-th block, I = yKJkp 
is the dimensionless stiffness parameter, v = v' /{L*u)) is 
the dimensionless loading rate, and 4'i = ^i/^o is the 
dimensionless friction force at the i-th block. 

The corresponding equation of motion of the 2D BK 
model is given in the dimensionless form by 



(43) 



where Ui^j = Uij/L* is the dimensionless displacement of 
the block {i,j). It is assumed here that the displacement 
of each block occurs only along the direction of the plate 
drive. The motion perpendicular to the plate motion is 
neglected. 

Often (but not always), the motion in the direction 
opposite to the plate drive is also inhibited by imposing 
an infinitely large friction for iii < (or Uij < 0) in 
either case of ID or 2D. It is also often assumed both 
in ID and 2D that the loading rate v is infinitesimally 
small, and put i/ = during an earthquake event, a very 
good approximation for real faults (Carlson et al., 1991). 
Taking this limit ensures that the interval time during 



successive earthquake events can be measured in units 
of irrespective of particular values of v. Taking the 
1/ — > limit also ensures that, during an ongoing event, no 
other event takes place at a distant place independently 
of this ongoing event. 



2. The friction law 

The friction force $ causing a frictional instability is a 
crucially important element of the model. Here, we refer 
to the following two forms for [A] a velocity weakening 
friction force (Carlson and Langcr, 1989a), and [B] a rate- 
and-state dependent friction (RSF) law (Dieterich, 1979; 
Ruina, 1983; Marone, 1998; Scholz, 1998; Scholz, 2002). 

[A] In this velocity-weakening friction force, one simply 
assumes that the friction force (j) = is a unique 

function of the block velocity iii. In order for the model 
to exhibit a frictional instability corresponding to earth- 
quakes, one needs to assume a velocity- weakening force, 
i.e., (j){ui) needs to be a decreasing function of iii. The 
detailed form of (j){ui) would be irrelevant. The form 
originally introduced by Carlson and Langer has widely 
been used in many subsequent works, that is (Carlson et 
al, 1991), 



4>Wi) = \ 1-5 

1+2q!M,/(1-(5) ' 



for Ui < 0, 
for iii > 0, 



(44) 



where the maximum value corresponding to the static 
friction has been normalized to unity. This normaliza- 
tion condition (/)(ui = 0) = 1 has been utilized to set 
the length unit L* . The friction force is characterized by 
the two parameters, S and a. The former, d, introduced 
in (Carlson et al.,1991) as a technical device facilitating 
the numerics of simulations, represents an instantaneous 
drop of the friction force at the onset of the slip, while the 
latter, a, represents the rate of the friction force getting 
weaker on increasing the sliding velocity. As emphasized 
by Rice (Rice, 1993), this purely velocity- weakening fric- 
tion law applied to the discrete BK model did not yield 
a sensible continuum limit. To achieve the sensible con- 
tinuum limit, one then needs to introduce an appropri- 
ate short-length cutoff by introducing, e.g., the viscosity 
term as was done in (Mayers and Langer, 1993): See also 
the discussion below in subsecton 6. 

We note that, in several simulations on the BK model, 
the slip-weakening friction force (Ida, 1972; Shaw, 1995; 
Myers et al., 1996), where the friction force is assumed 
to be a unique function of the slip distance 4){ui), was 
utilized instead of the velocity-weakening friction force. 
Statistical properties of the corresponding BK model, 
however, seem not so different from those of the velocity- 
weakening friction force. 

Real constitutive relations is of course more complex, 
neither purely velocity-weakening nor slip- weakening. As 
discussed in section II, the RSF friction law was intro- 
duced to account for such experimental features, which 
we now refer to. 
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[B] From Eq. (IT6l) , friction force in the BK model is given 
by 

4 = {c' + a'log(4)+&'log4?}-^. (45) 

where M is an effective normal load. See section II. C for 
the other quantities and parameters. Among the several 
evolution laws, we use the aging (slowness) law fEg.lfTT])'). 



dt' 



= 1 - 



C 



(46) 



Under the evolution law above, the state variable ■ grows 
linearly with time at a complete halt v[ = reaching a 
very large value just before the seismic rupture, while it 
decays very rapidly during the seismic rupture. 

The equation of motion can be made dimensionless 
by taking the length unit to be the characteristic slip 
distance £ and the time unit to be the rise time of an 
earthquake uj^^ = [ni/kpY/'^ . Then, one has. 



= {vt - Ui) + P {ui+i - 2ui + Ui-i) 

- {c + a\og{v,/v*) + b\og{v*e,)) (47) 

= l-vA (48) 



df^ 
dt 



where the dimensionless variables are defined hy t ~ ujt' , 
u, = u'JC, = v'JiCuj), V* = </(£w), 6, = uO'^, 
V = v'/{Cuj), a = a'Af/{kpC), b = b'Af/{kpC), c = 
c'J\f/ (kpC), while I = (kc/kp)^/^ is the dimensionless stiff- 
ness parameter defined above. In some numerical simu- 
lations, a slightly different form is used for the a-term, 
where the factor inside the a-term, v/v*, is replaced by 
1 -I- {v/v*), i.e., 



2ui + Ui-i) 



{c + alog{l + ^)+b\oge^), (49) 

V* 



where the constant factor c in Eq. (|49| is shifted by 
blogv* from c in Ea.(|I5|. 

This replacement enables one to describe the system at 
a complete halt, whereas, without this replacement, the 
system cannot stop because of the logarithmic anomaly 
occurring at w = 0. Similar replacement is sometimes 
made also for the &-term, i.e., 6 to 1-1-6'. 

The values of various parameters of the model describ- 
ing natural faults were estimated (Ohmura and Kawa- 
mura, 2007). Typically, u!~^ corresponds to a rise time 
of an earthquake event and is estimated to be a few sec- 
onds from observations. Though the characteristic slip 
distance C remains to be largely ambiguous, an estimate 
of order a few mm or cm was given by Tse and Rice 
(Tse and Rice, 1986) and by Scholz (Scholz, 2002). The 
loading rate associated with the plate motion is typi- 
cally a few cm/year, and the dimensionless loading rate 
1/ = v' ji^Lbj) is of order v ~ 10~*. The dimensionless 



quantity kpC/N was roughly estimated to be of order 
lO"**. The dimensionless parameter c should be of order 
10'^ ^ 10^, and the a and b parameters are one or two 
orders of magnitude smaller than c. 



3. The ID BK model with short-range interaction 

The simplest version of the BK model might be the ID 
model with only nearest-neighbor inter-block interaction. 
Since this model was reviewed in an earlier RMP review 
article by Carlson, Langer and Shaw, 1994, we keep the 
discussion here to be minimum, focusing mainly on recent 
results obtained after the above review article. 

Earlier studies on the ID BK model have revealed that, 
while smaller events persistently obeyed the GR law, i.e., 
staying critical or near-critical, larger events exhibited a 
significant deviation from the GR law, being off-critical 
or "characteristic" (Carlson and Langer, 1989a; Carlson 
and Langer, 1989b; Carlson et al., 1991; Carlson, 1991a; 
Carlson, 1991b; Schmittbuhl, Vilotte and Roux, 1996). 

In Fig.[9l we show the recent data of the magnitude dis- 
tribution (Mori and Kawamura, 2005; 2006). The mag- 
nitude of an event, M , is defined by 



M = In I ^ At 



(50) 



where the sum is taken over all blocks involved in the 
event. 

As can be seen from Fig. |9l the form of the calculated 
magnitude distribution R{M) depends on the value of 
the velocity-weakening parameter a. The data for a = 1 
lie on a straight line fairly well, apparently satisfying the 
GR law, which may be called "near-critical" behavior. 
The values of the exponent B describing the power-law 
behavior is estimated to be S ~ 0.50 corresponding to 
the 6-value, b = |i? ~ 0.75. By contrast, the data 
for larger a deviate from the GR law at larger magni- 
tudes, exhibiting a pronounced peak structure, while the 
power-law feature still remains for smaller magnitudes: 
See Fig. (HUa). These features of the magnitude distri- 
bution were observed in many simulations in common 
(Carlson and Langer, 1989a; Carlson and Langer, 1989b; 
Carlson et al, 1991). It means that, while smaller events 
exhibit self-similar critical properties, larger events tend 
to exhibit off-critical or characteristic properties, which 
may be called "supercritical" . Te data for smaller a < 1 
exhibit still considerably different behaviors from those 
for a > 1. Large events are rapidly suppressed, which 
may be called "subcritical" behavior. For a = 0.25, in 
particular, all events consist almost exclusively of small 
events only: See Fig. [5l^b). Here the words "critical", 
"supercritical" and "subcritical" have been defined on 
the basis of the shape of the magnitude-frequency rela- 
tionship. 

As an example of properties other than the magnitude 
distribution, we show in Fig. [lU] the recurrence-time dis- 
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FIG. 9: The magnitude distribution of earthquake 
events of the ID BK model with nearest-neighbor 
interaction for various values of the friction parameter 
a; (a) for larger a = 1, 2, 3, 5 and 10, and (b) for smaller 
a = 0.25,0.5,0.75 and 1. The parameters I and 5 are 
fixed to be Z = 3 and 5 = 0.01. The system size is 
N = 800. Taken from (Mori and Kawamura, 2006). 



tribution (Mori and Kawamura, 2005; 2006). The recur- 
rence time T is defined here locally for large earthquakes 
with M > Mc = 3 or Ale = 4, i.e., the subsequent large 
event is counted when a large event occurs with its epi- 
center in the region within 30 blocks from the epicenter 
of the previous large event. As can be seen from the fig- 
ure, the tail of the distribution is exponential at longer 
T irrespective of the value of a. Such an exponential tail 
of the distribution has also been reported for real seis- 
micity (Corral, 2004). By contrast, the distribution at 
shorter T is non-exponential and largely differs between 
for a = 1 and for a > 1. For a > 1, the distribution 
has an eminent peak corresponding to a characteristic 
recurrence time, which suggests the near-periodic recur- 
rence of large events. Such a near-periodic recurrence of 
large events was reported for several real faults (Nishenko 
and Buland, 1987; Scholz, 2002). For a =_1, by con- 
trast, the peak located close to the mean T is hardly 



discernible. Instead, the distribution has a pronounced 
peak at a shorter time, just after the previous large event. 
In other words, large events for a ~ 1 tend to occur as 
"twins" . A large event for the case of a = 1 often occurs 
as a "unilateral earthquake" where the rupture propa- 
gates only in one direction, hardly propagating in the 
opposite direction. 




T/T 

FIG. 10: The local recurrence-time distribution of the 
ID BK model with nearest-neighbor interaction for 
various values of the frictional parameter a. Large 

events of M > Mc = 3 or 4 are considered. The 
parameters are I and S are I = 3 and S — 0.01. The 
recurrence time T is normalized by its mean T. The 
total number of blocks is = 800. The insets represent 
the semi-logarithmic plots including the tail part of the 
distribution. Taken from (Mori and Kawamura, 2005). 

Possible precursory phenomenon exhibited by the 
model is of much interest, since it might open a way to an 
earthquake forecast. In fact, certain precursory features 
were observed in the ID BK model. Shaw, Carlson and 
Langer examined the spatio-temporal patterns of seismic 
events preceding large events, observing that the seismic 
activity accelerates as the large event approaches (Shaw, 
Carlson and Langer, 1992). Mori and Kawamura ob- 
served that the frequency of smaller events was gradually 
enhanced preceding the mainshock, whereas, just before 
the mainshock, it is suppressed in a close vicinity of the 
epicenter of the upcoming event (Mori and Kawamura, 
2005; 2006), a phenomenon closely resembling the "Mogi 
doughnut" (Mogi, 1969; 1979; Scholz, 2002). Fig.Ejrep- 
rcsents the space-time correlation function between the 
large events and the preceding events of arbitrary size 
(dominated in number by smaller events): It represents 
the conditional probability that, provided that a large 
event of M > Mc = 3 occurs at a time and at a spatial 
point ro , an event of arbitrary size occurs at a time to — t 
and at a spatial point Tq ± r. As can be seen from the 
inset of Fig. [TTJ seismic activity is gradually accerelated 
toward the mainshock either spatially or temporally. As 
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can be seen from the main panel, however, the seismic 
activity is supressed just before the mainshock in a close 
vicinity of the epicenter of the mainshock: See the dip 
developing around r = for t < 0.01. 

It turned out that the size of the quiescence region was 
always of only a few blocks, independent of the size of the 
upcoming mainshock (Mori and Kawamura, 2006). This 
may suggest that the quiescence is closely related to the 
discrete nature of the BK model: See subsection III. A. 6 
below. As such, the size of the quiescence region cannot 
be used in predicting the size of the upcoming mainshock. 
Instead, certain correlation was observed between the size 
of the upcoming mainshock and the size of the seismically 
active "ring" region surrounding the quiescence region 
(Pepke, Carlson and Shaw, 1994; Mori and Kawamura, 
2006). Such a correlation was also reported in real seismic 
catalog (Kossobokov and Carlson, 1995). 



;v=0-0.033 - 
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FIG. 11: The event frequency preceding the large event 
of M > Mc = 3 versus the distance from the epicenter 
of the upcoming mainshock of the ID BK model with 
nearest-neighbor interaction. The parameters a, / and d 
are a = 1, / = 3 and 5 ~ 0.01. The data are shown for 
several time periods before the mainshock. The insets 
represent similar plots with longer time intervals. The 
system size is iV = 800. Taken from (Mori and 
Kawamura, 2006). 

An aftershock sequence obeying the Omori law, al- 
though a common observation in real seismicity, is not 
observed in the BK model, at least in its simplest version 
(Carlson and Langer. 1989a, 1989b; Mori and Kawamura, 
2006). Interestingly, Pelletier reported that the inclusion 
of the viscosity effect in the form of "dashpot" in the 
2D BK model, together with the introduction of inhomo- 
geneity of friction parameters, could realize an aftershock 
sequence obeying the Omori law (Pelletier, 2000). The 
frictional force employed by Pelletier was a very simple 
one, i.e., a constant dynamical vs. static friction coeffi- 
cient. Further analysis will be desirable to establish the 
occurrence of the aftershock sequence obeying the Omori 
law in the BK model. 



We note in passing that the ID BK model has also 
been extended in several ways, e.g., taking account of 
the effect of the viscosity (Myers and Langer, 1993; Shaw, 
1994; De and Ananthakrisna, 2004; Mori and Kawamura, 
2008b), modifying the form of the friction force (Myers 
and Langer. 1993; Shaw, 1995; Cartwright, 1997; De and 
Ananthakrisna, 2004), and driving the system only at 
one end of the system (Vieira, 1992; 1996a). The effect 
of the long-range interactions introduced between blocks 
was also analyzed, which we will review in subsection 
III.A.4. 



4. The 2D BK model with short-range interaction 

Real earthquake faults are 2D rather than ID. Hence, 
it is clearly desirable to study the 2D version of the BK 
model in order to further clarify the statistical properties 
of earthquakes. The 2D BK model taken up here is to be 
understood as representing a 2D fault plane itself, where 
the direction orthogonal to the fault plane is not consid- 
ered explicitly in the model (Carlson, 1991b). The other 
possible version is the one where the second direction of 
the model is taken to be orthogonal to the fault plane 
(Myers et al, 1996). 

Extensive numerical studies have revealed that statis- 
tical properties of the 2D BK model are more or less 
similar to those of the ID BK model reviewed in the pre- 
vious subsection, at least qualitatively. The magnitude 
distribution R{M) of the 2D BK model was studied by 
several groups (Carlson and Langer, 1989a; Carlson and 
Langer, 1989b; Carlson ct al, 1991; Kumagai, et al, 1999; 
Mori and Kawamura, 2007). In Fig. [T^l we show typical 
behaviors of the magnitude distribution of the 2D BK 
model with varying the frictional parameter a (Mori and 
Kawamura, 2007). For smaller a ^ 0.5, R{M) bends 
down rapidly at larger magnitudes, exhibiting a "sub- 
critical" behavior. Only small events of M ^ 2 occur in 
this case. At a ^ 0.5, large earthquakes of their mag- 
nitudes M ~ 8 suddenly appear, while earthquakes of 
intermediate magnitudes, say, 2 ^ M ^ 6, remain rather 
scarce. Such a sudden appearance of large earthquakes at 
a = ad — 0.5 coexisting with smaller ones has a feature 
of a discontinuous or "first-order" transition. 

In this context, it might be interesting to point out that 
Vasconcelos observed that a single block system exhibited 
a "first-order transition" at a = 0.5 from a stick-slip to 
a creep (Vasconcelos, 1996), whereas this discontinuous 
transition becomes apparently continuous in many-block 
system (Vieira et al, 1993; Clancy and Corcoran, 2005). 
A "first-order" transition observed at a = Ud ~ 0.5 in 
the 2D model may have some relevance to the first-order 
transition of a single-block system observed by Vascon- 
celos, although events observed at a < ad in the present 
2D model are not really creeps, but rather are stick-slip 
events of small sizes. 

With increasing a further, earthquakes of intermediate 
magnitudes gradually increase their frequency. Fig. ll2f b) 



21 



exhibits R{M) for larger a. In the range of 1 5 a 5 10, 
R[M) exhibits a pronounced peak structure at a larger 
magnitude, deviating from the GR law, while it exhibits a 
near straight-line behavior corresponding to the GR law 
at smaller magnitudes ("supercritical" behavior). As a 
increases further, the peak at a larger magnitude becomes 
less pronounced. At a = ac2 — 13, R{M) exhibits a near 
straight-line behavior for a rather wide magnitude range, 
though R{M) falls off rapidly at still larger magnitudes 
M ^ 7, indicating that the "near-critical" behavior ob- 
served for a = ac2 — 13 cannot be regarded as a truly 
asymptotic one, since this rapid fall-off of R{M) at very 
large magnitudes is a bulk property, not a finite-size ef- 
fect. 



A "phase diagram" of the model in the elasticity pa- 
rameter I versus the friction parameter a, as reported by 
Mori and Kawamura, 2007 is shown in Fig. [T31 The re- 
gion or the "phase" , called "supercritical" , "near-critical" 
and "subcritical" are observed. The straight-line be- 
havior of R{M), i.e., the GR law is realized only in 
the restricted region in the phase diagram along the 
phase boundary between the supercritical and subcritical 
regimes. Even along the phase boundary, the GR relation 
is characterized by a finite cutoff magnitude above which 
larger earthquakes cease to occur. Hence, the GR rela- 
tion, as observed in a ubiquitous manner in real faults, 
is not realized in this model. Since each phase boundary 
has a finite slope in the a — I plane, one can also induce 
the "subcritical" -"supercritical" transition with varying 
the l-value for a fixed a (Espanol, 1994; Vieira, 1996b). 
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FIG. 12: The magnitude distribution R{M) of the 2D 
BK model with nearest-neighbor interaction for various 
values of the friction parameter a. The other 
parameters are / = 3 and 6 ~ 0.01. Fig. (a) represents 
R{M) for smaller values of the friction parameter 
< a < 3, while Fig.(b) represents R{M) for larger 
values of the friction parameter 3 < a < oo. The system 
size is 60 x 60. Taken from (Mori and Kawamura, 
2008a). 
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FIG. 13: Phase diagram of the 2D BK model with 
nearest-neighbor interaction in the friction parameter a 
versus the elastic-parameter I plane. The parameter S is 

6 = 0.01. Taken from (Mori and Kawamura, 2008a). 

As for other quantities, the recurrence-time distribu- 
tion of the 2D model exhibits a behavior similar to that 
of the ID model (Mori and Kawamura, 2007). As in case 
of ID, an aftershock sequence obeying the Omori law is 
not observed even in the 2D model, at least in its simplest 
version. The 2D model also exhibits precursory phenom- 
ena similar to the ones observed in the ID model (Mori 
and Kawamura, 2007). Acceleration of seismic activity 
prior to mainshock is observed in the supercritical regime, 
while it is not realized in the subcritical regimes. As in 
case of ID, mainshocks are accompanied by the "Mogi 
doughnut" -like quiescence in both supercritical and sub- 
critical regimes. 

As an other signature of the precursory phenomena, 
we show in Fig. [T3] the "time-resolved" local magnitude 
distribution calculated for time periods before the large 
event in the supercritical regime of a = 1 and / — 3 (Mori 
and Kawamura, 2007). Only events with their epicenters 
lying within 5 blocks from the upcoming mainshock of 
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its magnitude M > Mc = 5. As can be seen from the 
figure, an apparent -B- value describing the smaher magni- 
tude region gets smaller as the mainshock is approached, 
i.e., it changes from B ~ 0.89 of the long-time value to 
B ~ 0.65 in the time range ti^ < 0.1 before the main- 
shock. In real seismicity, an appreciable decrease of the 
_B-valuc has been reported preceding large earthquakes 
(Suyehiro, Asada and Ohtakc, 1964; Jaume and Sykes, 
1999; Kawamura, 2006). Obviously, a possible change 
in the magnitude distribution preceding the mainshock 
possesses a potential importance in earthquake fo recast. 




FIG. 14: The local magnitude distribution preceding 
the mainshock of M > Mc = 5 of the 2D BK model with 
nearest- neighbor interaction. The parameters are a = 1, 
I = 3 and 6 = 0.01. The data are shown for several time 
periods before the mainshock. The system size is 
60 X 60. Taken from (Mori and Kawamura, 2008a). 



5. The BK model with long-range interaction 

So far, we assumed that the interaction between blocks 
worked only between nearest-neighboring blocks. This 
may correspond to the situation where a thin isolated 
plate is subject to friction force and is driven by shear 
force (Clancy and Corcoran, 2006). However, a real fault 
is not necessarily a thin isolated plate, and the elastic 
body extends in a direction away from the fault plane. 
Indeed, the BK model extended in the direction orthogo- 
nal to the fault plane was also studied (Myers, Shaw and 
Langer, 1996). 

Considering the effect of such an extended elastic 
body adjacent to the fault plane under certain condi- 
tions amounts to considering the effective inter-block in- 
teraction to be long-ranged. Thus, taking account of the 
effect of long-range interaction might make the model 
more realistic. Rundle et al studied the properties of the 
2D cellular automaton version of the BK model with the 
long-range interaction decaying as 1/r^ (Rundle, et al. 



1995). Xia et al studied the ID BK model with a vari- 
able range interaction where a block is connected to its 
R neighbors with a rescaled spring constant proportional 
to 1/R (Xia et al, 2005; Xia et al, 2007). The type of 
the long-range model considered by Xia et al may be re- 
garded as a mean-field type, since the model reduces to 
the mean-field infinite-range model in the i? — > oo limit. 

One can also derive the relevant long-range interaction 
based on an elastic theory (Mori and Kawamura, 2008a) . 
Suppose that the 3D elastic body in which the 2D BK 
model lies is isotropic, homogeneous and infinite, and a 
fault surface is a plane lying in this elastic body and 
slips along one direction only. Then, a static approxima- 
tion for an elastic equation of motion for the elastic body 
would give rise to a spring constant between blocks decay- 
ing with their distance r as 1/r^. This static assumption 
is justified when the velocity of the seismic-wave propa- 
gation is high enough compared with the velocity of the 
seismic-rupture propagation. 

Properties of the 2D BK model with the long-range 
power-law interaction derived from an elastic theory, i.e., 
the one decaying as 1/r^, was investigated (Mori and 
Kawamura, 2008a). The interaction between the two 
blocks at sites (i,j) and {i',j') is given in the dimen- 
sionless form by 

which falls off with distance r as 1/r^. Then, the dimen- 
sionless equation of motion of the 2D long-range can be 
written as 

(52) 

If one restricts the range of interaction to nearest neigh- 
bors and takes the spatially anisotropic spring con- 
stant to be isotropic, Ix = Iz = I, one recovers the 
isotropic nearest-neighbor model described by Eq. 1421 
The "isotropy" assumption l^ = Iz is equivalent to 
putting the Lame's constant to vanish. In fact, in the 
short-range model, such a spatial anisotropy of the 2D 
BK model turned out to hardly affect the statistical prop- 
erties of the model in the sense that the properties of the 
anisotropic model was quite close to the corresponding 
isotropic model characterized by the mean spring con- 
stant I = {Ix + Iz) /'2 (Mori and Kawamura, 2008a). 

One might also consider the ID BK model with the 
long-range interaction (Mori and Kawamura, 2008a). 
One possible way to construct the ID model might be 
to impose the condition on the corresponding 2D model 
that the systems is completely rigid along the z-direction 
corresponding to the depth direction, i.e., u{x,z,t) = 
u{x,t). This yields an effective inter-block interaction 
decaying with distance r as l/r^. 
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with the dimcnsionlcss equation of motion 



(54) 



In Figs. [TST a) and (b), we show the magnitude distri- 
bution R{M) of the long-range 2D BK model for smaller 
and larger values of a, i.e., (a) < a < 10 and (b) 
10 < a < oo (Mori and Kawamura, 2008a). Similarly to 
the short-range case, three distinct regimes are observed 
depending on the a-value. The intermediate- a region 
corresponds to the supercritical regime where R{M) ex- 
hibits a pronounced peak at a larger magnitude, show- 
ing a characteristic behavior. Major difference from the 
short-range case is that the subcritical behavior realized 
in the short-range model in the smaller- and larger-a re- 
gion is now replaced by the near-critical behavior in the 
long-range model. Namely, for smaller a < ad ~ 2 and 
for larger a > ac2 ^ 25, R{M) exhibits a near straight- 
line behavior over a rather wide magnitude range, and 
drops off sharply at larger magnitudes. The associated 
S-value is estimated to be i? ~ 0.59 (a < ad) and 
B ~ 0.55 (a > ac2), which is rather insensitive to the 
a-value. This straight-line behavior of R{M) cannot be 
regarded as a truly critical one, since R[M) drops off 
sharply at very large magnitudes. As in the short-range 
case, the change from the supercritical to the near-critical 
behaviors at a = ac2 ~ 25 is continuous, while it is dis- 
continuous at a = ad — 2. 

Such a near-critical behavior realized over a wide pa- 
rameter range is in sharp contrast to the behavior of the 
corresponding short-range model where R{M) at smaller 
and larger a exhibits only a down-bending subcritical be- 
havior, while a straight-line near-critical behavior is re- 
alized only by fine-tuning the a-value to a special value 
a ~ ac2 ■ The robustness of the near-critical behavior of 
R{M) observed in the 2D long-range model might have 
an important relevance to real seismicity, since the GR 
law is ubiquitously observed for different types of faults. 
Note also that the associated B- value observed here turns 
out to be close to the one observed in real seismicity (Mori 
and Kawamura, 2008a). 

In Fig. 1161 the behavior of R{M) is summarized in 
the form of a "phase diagram" in the friction parameter 
a versus the elastic-parameter I plane (Mori and Kawa- 
mura, 2008a). As can be seen from the figure, the phase 
diagram of the long-range model consists of three distinct 
regimes, two of which are near-critical regimes and one is 
a supercritical regime. The "phase boundary" between 
the smaller- a near-critical regime and the supercritical 
regime represents a "discontinuous transition" , while the 
one between the larger-a near-critical regime and the su- 
percritical regime represents a "continuous transition" . 
For comparison, the corresponding phase boundary of 
the short-range model is also shown. The near-critical 
phases in the long-range model are replaced by the sub- 
critical phases in the short-range model. 

It might be interesting to notice that the system at 
different "phases" of Fig. [1^] really show different prop- 
erties. For example, we show in Fig. [T7] the magnitude 




FIG. 15: The magnitude distribution R{M) of the 2D 
BK model with long-range interaction for various values 
of the friction parameter a. The other parameters are 
/ — 3 and 5 = 0.01. Fig. (a) represents R{M) for smaller 
values of the frictional parameter < a < 10, while 

Fig.(b) represents R{M) for larger values of the 
frictional parameter 10 < a < oo. The system size is 
60 X 60. Taken from (Mori and Kawamura, 2008a). 



dependence of the mean displacement Au at a seismic 
event (Mori and Kawamura, 2008a) . As can be seen from 
the figure, the data in the two near-critical regimes (the 
data in blue and in green) are grouped into two distinct 
branches, while the data in the supercritical regime (the 
data in red) exhibit a significantly different behavior. In- 
terestingly, the mean displacement in the near-critical 
regimes hardly depends on the event magnitude. 

It was observed that the mean stress drop at a seis- 
mic event also hardly depends on the event magnitude in 
the near-critical regimes of the 2D long-range BK model 
(Mori and Kawamura, 2008a). A similar independence 
was also reported in the mean-field-type ID long-range 
BK model (Xia et al, 2005; 2008) and in the ID long- 
range BK model (Mori and Kawamura, 2008a). 
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"continuous 
Transition" 




FIG. 16: The phase diagram of the 2D BK models with 
long-range interaction in the friction parameter a versus 
elastic-parameter I plane, which is compared with the 
one of the 2D BK model with short-range interaction. 
The parameter 8 is set 5 ~ 0.01. Taken from (Mori and 
Kawamura. 2008a). 
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FIG. 17: The magnitude dependence of the mean 
displacement Au at each seismic event of the 2D BK 
model with long-range interaction. In the main panel, 
the friction parameter a is varied with fixing the system 
size 60 X 60, while in the inset the system-size iV is 
varied for the case of a = 30. The parameters / and 8 
are fixed tol ~?> and 5 ~ 0.01. Taken from (Mori and 
Kawamura, 2008a). 



oj] scale into the problem. Therefore, in order to check 
the validity of the model, it is crucially important to ex- 
amine the continuum limit of the BK model carefully. 
Indeed, Rice criticized that the discrete BK model with 
the velocity-weakening friction law was "intrinsically dis- 
crete" , lacking in a well-defined continuum limit (Rice, 
1993). Rice argued that the spatiotemporal complex- 
ity observed in the discrete BK model was due to the 
inherent discreteness of the model, which should disap- 
pear in continuum. Indeed, he applied the RSF law, 
which possessed an intrinsic length scale corresponding 
to the characteristic slip distance, and showed that the 
system tended to exhibit a quasi-periodic behavior, if the 
grid spacing d' was taken smaller than the characteristic 
length scale, while if the grid spacing d' was taken longer 
than it, the system exhibited an apparently complex or 
critical behavior. This problem of the continuum limit of 
the BK model was also addressed by Myers and Langer 
(Myers and Langer, 1993) within the velocity-weakening 
friction law, who introduced the Kelvin viscosity term 
to produce a small length scale which allowed a well- 
defined continuum limit. Myers and Langer, and sub- 
sequently Shaw (Shaw, 1994), observed that the added 
viscosity term smoothed the rupture dynamics, appar- 
ently giving rise to the continuum limit accompanied by 
the spatiotemporal complexity. More recently, the con- 
tinuum limit of the ID BK model with and without the 
viscosity was examined by Mori and Kawamura within 
the velocity- weakening friction law (Mori and Kawamura, 
2008b). 

Thus, two different ways of taking the continuum limit 
of the BK model were tried so far, each introducing the 
short length scale via (A) the viscosity term, or (B) the 
RSF law. In this subsection, we examine the former (A), 
while the latter (B) will be discussed in the next subsec- 
tion. 

As mentioned, the naive continuum limit of the dis- 
crete BK model with a velocity-weakening friction force 
without viscosity has a problem in that the pulse of slip 
tends to become increasingly narrow in width in the limit, 
i.e., the dynamics becomes sensitive to the grid spacing 
d' — > 0. One way to circumvent this problem is to in- 
troduce the viscosity term ri'd^Ui/{dx'^dt') into EgHT] 
to produce a small length scale, where 77' is the viscosity 
coefficient. Myers and Langer showed that, owing to the 
added viscosity term, the system became independent of 
the grid spacing d' as long as a new small length scale e', 
defined by 



(55) 



6. Continuum limit of the BK model 

Although the BK model has widely been used as a useful 
tool to investigate statistical properties of earthquakes, 
the block discretization inherent to the model construc- 
tion is a crude approximation of the originally contin- 
uum earthquake fault. It introduces the short-length cut- 



is sufficiently larger than the grid spacing d' (Myers and 
Langer, 1993). With ^' being the wave velocity in the 
continuum limit, this small length scale e' can also be 
given in the dimensionless form as 



e EE £7(^/0.) = 7rJ^, 
V a 



(56) 
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where 77 = rj' / {^^'^ / u>) is the dimensionless viscosity coef- 
ficient. The dimensionless distance r between the block 
i and i' is measured by 

r = d\i-i'l (57) 

where d = d' /lj) is the dimensionless grid spacing. 
The continuum limit corresponds to taking the limit 
d — >■ with fixing L = Nd and r, which means N 00 
and I — >■ 00. Thus, taking the continuum limit in the 
BK model corresponds to making the model to be in- 
finitely rigid / — (X). Numerically, various observables 
were calculated with successively smaller d to examine 
its asymptotic c? — > limit. 

Shaw showed, by adding the viscosity term to the ID 
BK model, that the magnitude distribution became in- 
dependent of the grid spacing d' for sufficiently small 
d' (Shaw, 1994). Mori and Kawamura studied the ID 
BK model with successively smaller grid spacings d' to 
examine how various statistical properties of the model 
changed and approached the continuum limit for both 
cases of nonzero (77 > 0) and zero (77 = 0) viscosity (Mori 
and Kawamura, 2008b). It was then observed that, in the 
former viscous case, the results converged to the contin- 
uum limit when the condition d < e was met, whereas, 
in the latter non-viscous case, such a convergence was 
obscure. 

As an example, we show in Fig. [15] the way of con- 
vergence of the magnitude distribution function R{M) 
for a = 1 (a) and for a = 3 (b), in the viscous case 
(77 = 0.02). For both cases of a = 1 and 3, the continuum 
limit seems to be well reached, i.e., R{M) seems to con- 
verge to an asymptotic form for smaller d, except that the 
minimum magnitude continuously gets lower as the grid 
spacing d gets smaller. A similar result was reported by 
Shaw, 1994. From Fig.fTSla'l. one also sees that a nonzero 
viscosity tends to weaken the GR character of the mag- 
nitude distribution somewhat. Such a deviation from 
the GR law at smaller magnitudes is probably originated 
from the fact that the viscosity tends to make the rela- 
tive displacement of neighboring blocks being smoother, 
enhancing the correlated motion of neighboring blocks, 
which makes the frequency of smaller events of one or a 
few blocks considerably reduced (Mori and Kawamura, 
2008b). 

The small-length cutoff scale e as given by Eq. [55] is 
estimated here to be e ~ 0.44 and 0.26 for a = 1 and 
3, respectively. As can be seen from Figs. [TSTa) and (b), 
R{M) converges to an asymptotic form for the a-values 
smaller than d ~ 1/4 and 1/8 for a = 1 and 3, respec- 
tively, which is consistent with the expected condition of 
the continuum limit d < e. 

As mentioned in subsection III. A. 3, the BK model gen- 
erally gives rise to a seismic quiescence phenomenon prior 
to mainshock, i.e., the Mogi-doughnut. Then, a natu- 
ral question is whether the doughnut-like quiescence ob- 
served in the discrete BK model survives the continuum 
limit, or it is a phenomenon intrinsically originated from 
the short cutoff length scale of the model. This question 




FIG. 18: The magnitude distribution R{M) of 
earthquake events of the ID viscous BK model 
(77 = 0.02) with 5 = 0.01. The dimensionless grid 
spacing d is varied in the range 1 > d > 1/32. Figs. (a) 
and (b) represent the cases of a = 1 and 3, respectively. 
The system size is L = dN = 200. Taken from (Mori 
and Kawamura, 2008b). 



was addressed in (Mori and Kawamura, 2008b). Fig. [T9l 
exhibits the time-dependent spatial correlation functions 
before the mainshock in the case of the viscous model of 
a — 1. As the grid spacing d gets smaller, the spatial 
range of the quiescence gets narrower, tending to vanish 
for small enough d: See the inset of Fig. [TO] This observa- 
tion strongly suggests that the doughnut-like quiescence 
might vanish altogether in the continuum limit d — ^ 0. 
Thus, the doughnut-like quiescence observed in the dis- 
crete BK model is likely to be a phenomenon closely re- 
lated to the short-length cutoff scale of the model. This 
seems fully consistent with the observation that the one- 
block events are responsible for the observed doughnut- 
like quiescence (Mori and Kawamura, 2006; 2008a). 

The observation might have some implications to real 
seismicity. While the real crust is obviously a contin- 
uum, it is often not so uniform, possibly with a short- 
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FIG. 19: The event frequency in the time period 
tv ~ Q ^ 0.01 immediately before the mainshock of 
M > Mc = 2 of the ID viscous BK model (r/ = 0.02) 
with a = 1 plotted versus r, the distance r from the 
epicenter of the upcoming mainshock. The 
dimensionless grid spacing d is varied in the range 
1/4 > d > 1/32. The parameter 5 is fixed to 5 = 0.01. 
The system size is L = dN = 200. The insets represent 
the peak position of the event frequency, corresponding 
to the range of the doughnut-like quiescence, as a 
function of the dimensionless grid spacing d. The 
doughnut-like quiescence vanishes in the continuum 
limit d — > 0. Taken from (Mori and Kawamura, 2008b). 



length cutoff. In any case, in real earthquakes the Mogi- 
doughnut is occasionally reported to occur (Mogi, 1969; 
1979; Scholz, 2002), although establishing its statistical 
significance is sometimes not easy. Then, our present 
result may suggest that, if the real crust possesses a cut- 
off length scale due to the inhomogeneity of the crust, 
the "Mogi-doughnut" quiescence might occur at such a 
length scale. In other words, spatial inhomogeneity might 
be an essential ingredient for the Mogi-doughnut to occur 
in real seismicity (Mori and Kawamura, 2008b). 



7. The BK model with RSF law 

So far, wc have mostly assumed a simple velocity- 
weakening friction law where the friction force is a single- 
valued function of the velocity. As detailed in section II 
and in subsection III. A. 2, the RSF law is now regarded 
in seismology as the standard consititutive law. 

Tse and Rice employed this RSF constitutive relation 
in their numerical simulations of earthquakes (Tse and 
Rice, 1986). These authors studied the stick-slip motion 
of the two-dimensional strike-slip fault within an elas- 
tic continuum theory, assuming that the fault motion is 
rigid along strike. It was then observed that large events 
repeated periodically. Since then, similar RSF consti- 



tutive laws have widely been used in numerical simu- 
lations (Stuart, 1988; Horowitz and Ruina, 1989; Rice, 
1993; Ben-Zion and Rice, 1997; Kato and Hirasawa, 1999; 
Kato, 2004; Bizzarri and Cocco, 2006). Somewhat dif- 
ferent type of slip- and state-dependent constitutive law 
was also used (Cochard and Madariaga, 1996). 

Cao and Aki performed a numerical simulation of 
earthquakes by combining the ID BK model with the 
RSF law in which various constitutive parameters were 
set nonuniform over blocks (Cao and Aki, 1986). Ohmura 
and Kawamura extended an earlier calculation by Cao 
and Aki to study the statistical properties of the ID BK 
model combined with the RSF constitutive law with uni- 
form constitutive parameters (Ohmura and Kawamura, 
2007). Clancy and Corcoran also performed a numeri- 
cal simulation of the ID BK model based on a modified 
version of the RSF law (Clancy and Corcoran, 2009). 

Rice and collaborators argued that the slip complexity 
of the BK model might be caused by its intrinsic discrete- 
ness (Rice, 1993; Ben-Zion and Rice, 1997). In this con- 
text, it is important to clarify the statistical properties of 
the model where the discrete BK structure is combined 
with the RSF law, to compare its statistical properties 
with those of the standard BK model with the velocity- 
weakening or slip- weakening friction law reviewed in the 
previous subsections. 

Recent study by Morimoto and Kawamura has re- 
vealed that the model exhibits largely different behaviors 
depending on whether the frictional instability is either 
"strong" or "weak" (Morimoto and Kawamura, 2011). 
The condition of strong or weak frictional instability is 
given by 6 > 2/^-1-1 or 6 < 2^^-1-1, respectively, for the ID 
BK model. In the case of a weaker frictional instability, 
the model exhibits a precursory process where a slow nu- 
cleation process occurs prior to mainshock. In the next 
subsection, we discuss such a precursory process realized 
in the BK model in more detail. Interestingly, presence 
or absence of such a nuclcation process also affects sta- 
tistical properties of the model. From a simulation point 
of view, the case of a weaker friction instability is much 
harder to deal with, since slow and long-standing nucle- 
ation process prior to mainshock generally requires a lot 
of CPU time. 

Statistical properties of the ID BK model with the 
RSF law EqEl(or Eqgni) and Eq|15]was investigated by 
Ohmura and Kawamura for the case of a strong frictional 
instability (Ohmura and Kawamura, 2007), and by Ya- 
mamoto and Kawamura for the case of a weak frictional 
instability (Yamamoto and Kawamura, 2011). Typical 
behaviors of the magnitude distribution are respectively 
shown in Figs. [20la) and (b). As can be seen from the 
figure, when the frictional instability is strong, almost 
flat distribution spanning from small to large magnitudes 
is realized, while, as the critical value is approached, a 
peak at a larger magnitude becomes more pronounced 
giving rise to an enhanced characteristic behavior. In the 
weak frictional instability regime, the distribution has no 
weight at smaller magnitudes, with a pronounced peak 
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only at a large magnitude. It means that only large earth- 
quakes of more or less similar magnitude occ ur in the 
regime of a weak frictional instability. 




I Z 3 4 



M 

FIG. 20: (Color online) The magnitude distribution of 
the ID BK model with the RSF law, for the case of (a) 
a strong frictional instability b > be, and of (b) a weak 

frictional instability b < be, with be = 2P + 1. The 
parameter values are a = 0, c = 1000, u = 10^^, v* = 1 
and / = 3 in (a), and a = 1, 6 = 5, c = 1000 v* = 1 and 
I = 5 in (b). The borderline 6- value is be = 19 in (a), 
and 6c = 51 in (b). The system size is = 800 in (a), 
and N = 1200 in (b). (a) Taken from (Ohmura and 
Kawamura, 2007). (b) Taken from (Morimoto and 
Kawamura, 2011). 

Statistical properties of the corresponding 2D model 
were investigated by Kakui and Kawamura for both cases 
of weak and strong frictional instabilities (Kakui and 
Kawamura, 2011). In the 2D BK model, the condition of 
strong or weak frictional instability is given by 6 > 4P + 1 
or b < 4/^ + 1, respectively. Typical behaviors of the mag- 
nitude distribution are shown in Figs. Ella) and (b) for 
the cases of strong and weak instabilities, respectively. 
As can be seen from the figure, when the frictional insta- 
bility is strong, a behavior more or less close to the GR 
law, characterized by the exponent close to B ^ 2/3, is 



realized, although there is a weak shoulder-like structure 
superimposed at larger magnitudes. The observation of 
a near-critical behavior close to the GR law would be of 
much interest in conjunction with real seismicity. As the 
critical value is approached, on the other hand, a peak at 
a larger magnitude is further developed, giving rise to an 
enhanced characteristic behavior. In the weak frictional 
instability regime, the distribution has double peaks ex- 
hibiting more characteristic behavior: See Fig. [2lT bV 
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FIG. 21: (Color online) The magnitude distribution of 
the 2D BK model with the RSF law, for the case of (a) 
a strong frictional instability b > be, and of (b) a weak 

frictional instability b < be, with be = + 1. The 
parameter values are a — 1, c — 1000, ly ~ 10^*, w* = 1 
and / = 2 in (a), and a = 1, c = 1000, v = 10"^, u* = 1 
and Z = 2 in (b). The borderline value is 6c = 17 in both 
(a) and (b). The system size is = 60 x 60 in (a), and 

= 30 X 30 in (b). Taken from (Kakui and Kawamura, 
2011). 



8. Nucleation process of the BK model 

In this subsection, we touch upon the nucleation process 
as a precursory phenomenon prior to mainshock as re- 
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alized in the BK model obeying the RSF law. It was 
observed that the nucleation process is realized even in 
the BK model with the RSF law for both cases of ID and 
2D, if the model lies in the regime of a weak frictional 
instability (Morimoto and Kawamura, 2011; Kakui and 
Kawamura, 2011). Namely, prior to seismic rupture, the 
system exhibits a slow rupture process localized to a com- 
pact "seed" area with its rupture velocity orders of mag- 
nitude slower than the seismic wave velocity. The system 
spends a very long time in this nucleation process, and 
then at some stage, exhibits a rapid acceleration process 
accompanied by a rapid growth of the rupture velocity 
and a rapid expansion of the rupture zone, finally getting 
into a final seismic rupture or a mainshock (Dieterich, 
2009). Such a nucleation process has also been observed 
and extensively studied in the continuum model: See, 
e.g., (Ampuero and Rubin, 2008). We illustrate in Fig. [22] 
typical example of seismic events realized in the ID BK 
model with the RSF law for each case of a weak frictional 
instability (b), and of a strong frictional instability (a). 
As can be seen from the figure, a slow nucleation process 
with a long duration time is observed only in (b), while 
such a nucleation process is absent in (a). 
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As mentioned, the condition for the appearance of such 
a nucleation process is given by 6 < &c = 2/^ -I- 1 in ID, 
and by 6 < 6c = 4^^ + 1 in 2D (for a square array of 
blocks). Indeed, Morimoto and Kawamura found that 
the critical nucleation size at which the slow nucleation 
process ends getting into the acceleration stage is given 
by Xc = 7r/[arccos(l — ^jr)] — 1 in units of block size 
(Morimoto and Kawamura, 2011). Indeed, this length ATc 
corresponds in its physical meaning to the length h* of 
Rice (Rice, 1993), although its detailed functional form, 
e.g., the dependence on 6, is somewhat different from the 
standard one. The condition of this critical nucleation 
size being greater than the block size Xc > 1 yields the 
condition of the weak frictional instability b < be. In 
other words, when b > b^ the nucleation process cannot 
be realized in the BK model due to its intrinsic discrete- 
ness. Indeed, this is exactly the situation as discussed by 
Rice (Rice, 1993). 

The above observation means that, if one takes the 
continuum limit of the BK model with the RSF law, 
the system should necessarily lie in the limit of a weak 
frictional instability, since the continuum limit means 
I ^ CO. Hence, at least as long as one considers a uni- 
form fault obeying the RSF law without any discretization 
short-length scale, earthquakes should exhibit character- 
istic properties rather than critical properties. This fully 
corroborates an earlier criticism by Rice against the SOC 
view of earthquakes based on the BK model (Rice, 1993). 
Indeed, in seismology the concept of earthquake cycle has 
been used in long-term probabilistic earthquake forecasts 
(Scholz, 2002; Nishenko, 1987; Working Group on Cali- 
fornia Earthquake Probabilities, 1995). Of course, a big 
issue to understand is what is then the true origin of the 
GR law widely observed in real seismicity. 




B. Continuum models 
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FIG. 22: (Color online) The typical rupture process 
realized in the ID BK model with the RSF law for (a) a 
strong and (b) a weak frictional instability, each 
corresponding to (a) b > be and (b) b < be with 
be = 2P -\- 1. The color represents the rupture velocity. 
The parameter values are a = 1, c = 1000, ly = 10~^, 
V* ~ 1 and / = 5 for both (a) and (b) corresponding to 
be = 51, whereas 6 = 60 in (a) and 6 = 3 in (b). Taken 
from (Morimoto and Kawamura, 2011). 



As discussed in III. A. 6. [Ricd ( 1993() criticized inherently 
discrete models, where simulated earthquake sequences 
depend on computation grid size. He confirmed in nu- 
merical simulations that complex earthquake sequences 
disappear when the grid size is sufficiently smaller than 
the critical size of slip nucleation zone for almost spa- 
tially uniform frictional properties. Moreover, he argued 
that geometrical and/or material disorder is the origin of 
complexity of earthquakes. The models with sufficiently 
small grid sizes may be called continuum models, which 
generate simulation results independent of the grid size, 
in contrast to inherently discrete models. Note that if a 
model does not have a finite critical size for nucleating 
unstable slip, such as a model with constant static and 
dynamic friction, it is always inherently discrete. In this 
subsection, we discuss continuum models of earthquakes, 
especially models using the rate- and state-dependent 
friction (RSF) law. In the RSF law, the critical size of 
slip nucleation can be defined as a function of frictional 
constitutive parameters, and the computation grid sizes 
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are sufficiently smaller than the critical size in the stud- 
ies mentioned below. We use elastic continuum models 
below, in contrast to spring-block models in the previous 
section. " Continuum model" is thus used to express two 
senses. 

The RSF law has commonly been used in mod- 
els for und erstanding earthquake phenomena ( Dieterichl . 
120091 : [Schok . ,20021) .These models were sometimes con- 
structed for reproducing and understanding particular 
earthquakes, earthquake cycles, or sliding processes ob- 
served by seismometers, strainmeters. Global Positioning 
System (GPS), etc. We will see deterministic aspects of 
earthquake phenomena, in addition to statistical char- 
acteristics of earthquake s. Note that comprehensive re- 
views were presented b y Ben-ZionI ( 2008() : iRundle et all 
(|2003l ): iTurcotte efall (|2009( ) for models of statistical 
properties of earthquakes using friction laws other than 
the RSF law. 



1. Earthquake cycles, asperities, and aseismic sliding 

Before introducing earthquake models, we briefly re- 
view observational facts about earthquakes and fault slip 
behavior. Earthquakes repeatedly occur at the same 
fault segment. At the Parkfield segment along the San 
Andreas fault, California, magnitude of about 6 inter- 
plate earthquakes have occurr ed at recurrence i n terval s 
of 23 ± 9 years since 1857 (jSvkes and Menkd . |2006[ ). 
Great earthquakes of magnitude 8 class repeatedly oc- 
curred along the Nankai trough, where the Philippine 
Sea plate subducts beneath southwest e rn Jap an, every 
one hundred years (|Svkes and Menkd . l2006li . Quasi- 
periodic earthquake recurrence has been used for long- 
term forecasts of earthquakes (Working Group on Cali- 
fornia Earthquake Probabilities, 1995; Matthews ct al., 
2002). One of the most remarkable examples of reg- 
ularity of earthquakes was found off Kamaishi, where 
the Pacific plate subducts beneath northern Honshu, 
Japan. Magnitude of 4.8 ± 0.1 earthquakes have repeat- 
edly occurred at recurrence inte rvals of 5.5 ± 0.7 y ears 
at the same region since 1957. lOkada et all ( 20031 ) es- 
timated coseismic slip distributions of recent Kamaishi 
earthquakes from seismic waveform data and found that 
they overlap with each other (Fig. [231) ■ Although 
many smaller earthquakes occur around the source area 
of the Kamaishi earthquakes, no comparable or larger 
earthquakes occur there. This observation suggests that 
aseismic sliding surrounds the source area of the Ka- 
maishi earthquakes, where stick-slip motion occurs, and 
steady loading by the surrounding aseismic sliding to 
the source area leads to the quasi-periodic recurrence 
of almost the same magnitude earthquakes. The vari- 
ance in the recurrence interval was suggested to come 
from temporal variation of as eismic sliding rate su rround- 
ing the earthquake source ( Uchida et all . l2005l ). Sig- 
nificant aftcrslip of the 2011 great Tohoku-oki earth- 
quake (M=9.0) rapidly loaded the source area of the 



Kamaishi earthquake, generating earthquakes at much 
shorter recurrence intervals. Recurrences of small earth- 
quakes at the same source areas in mainly creeping 
(aseissmic sliding) regions were found in many places and 
these earthquakes are called small repeating earthquakes 
(jlgarashi et all . l2003l : iNadeau and Johnson! . Il998t) . Al- 
though small earthquakes occur, most strain is released 
by aseismic sliding on these fault planes. The seismic 
coupling coefficient is defined by the long-term average 
of the ratio of seismic slip amount to total (seismic and 
aseismic) slip expected from relative plate motion. The 
seismic coupling coefficient is variable, dependent on lo- 
calities. It is close to unity at some segments along 
Chile and Aleutians, indicating little aseismic sliding and 
nearly complete locking during interseismic periods, and 
is nearly equal to zero at Maria nas, indicating no o r 
few large interplate earthquakes ( Pacheco et al . Il993f ). 



These facts show that aseismic sliding is common phe- 
nomenon and it plays an important part in strain release 
at plate boundaries and that frictional properties differ 
from place to place. 

A patch where stick-slip motion occurs, that is, a fault 
region where earthquakes repeatedly occur, is often called 
an asperity, which comes from the rock mechanics term 
for a contact spot between sliding surfaces as used in II. C. 
Note that an asperity of an earthquake occupies a con- 
siderable part of the earthquake fault area and its size is 
orders of magnitude larger than seismic slip amount. In 
contrast, an asperity of a sliding surface is much smaller 
and its size may be comparable to slip amount. The 
asperity model has been developed for explaining spa- 
tial heterogeneity in seismic sli p on faults and complex 
source processes of earthquakes ( Kanamori and McNallvl . 



19901) . When the asper- 



119821 : [Lav et al. l. ll982l : lThatcheif ] 

ity model was developed around 1980, sliding behavior 
surrounding asperities was not clarified from observations 
because aseismic sliding cannot be detected by seismome- 
ters. To detect aseismic sliding, geodetic observations 
such as GPS are require d. Since dense GPS ne tworks 
were established in 1990s (jSegaU and Davi£lll997^ ■ many 
aseismic sliding phenomena have been reported such as 
aftcrslip (postseismic sliding) and slow (silent) earth- 
quakes. The source areas of aftcrslip are usually located 
near coseismic slip areas (asperities), and the afterslip 
area and the asperity do not overlap a s shown in Fig. [Ml 



([Johnson et all [20061: [Mivazaki etdl . [2004 [Yagi et al 



[2003 ). which also support spatial heterogeneity of fric- 
tional properties. The locations of asperities of large 
earthquakes were confirmed to be locked during inter- 
seismic periods from geodetic observations ([Chlieh et all . 
[2008[ : [Hashimoto et all . [2009[ : [Perfettini et ad[2010f) . For 
instance, Figure [^ clearly shows that seismic slip areas 
of large interplate earthquakes off the Sumatra island co- 
incide with the locked areas during interseismic periods. 
For the 2011 great Tohoku-oki earthquake (M = 9.0), a 
significant peak of seismic slip larger than 30 m was esti- 
mated fi^HLiBX^ISio'^S-SLS^^^'^^'' waveform and tsunami 
data ( Koketsu et all . [201l[ ). This also suggests nonuni- 
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FIG. 23: (a) Recurrence of Kamaishi earthquakes of 
nearly the same magnitudes and recurrence intervals, 
(b) Cumulative seismic moment of Kamaishi 
earthquakes, (c) Coseismic slip distribution of the 1995 
and 2001 Kamaishi earthquakes estimated from seismic 
waveforms. Red broken contours and blue contours 
denote seismic slip of t he 1995 and 2001 e arthquakes, 
respectively ( Okada et aZ.I . |2Q03|) . 



form frictional property on the plate interface. 

Spatial distribution of asperities on plate boundaries 
has been estimated from source areas of past large inter- 
plate earthquakes, and e arthquakes repeatedly o c curre d 
on the same asperities ( Yamanaka and Kikuchil . 12004) . 
This suggest that the locations of asperities are un- 
changed at least a few earthquake cycles. Apparently 
complex earthquake cycle, where earthquake rupture ar- 
eas are variable, may be understood by a change in com- 
bination of simultaneously ruptured asperities. For ex- 
ample, two adjacent asperities are simultaneously rup- 
tured, resulting in a large earthquake in some cases, and 
one of them is ruptured to generate a smaller event in 
other cases. Note that some researchers object against 
persisten t asperities on the b asis of seismic waveform 
analyses (jPark and Moril2007l) . 
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FIG. 24: Spatial distribution of cumulative slip for 30 
days of aftcrslip of the 2003 Tokachi-oki earthquake (M 
= 8.0), off Hokkaido, northern Japan, estimated from 
GPS data (color contours) bv iMivazaki et all (|2004[) . 
Black contours with 0.5m interval show seismic slip in 
the 1973 Nemuro-oki (right), 1968 Tokachi-oki (left), 

and 2003 Tokachi-oki (cen ter) earthquakes 
( Yamanaka and Kikuchi |2004|) . The black star and 
small circles denote the epicenter and aftershocks of the 
2003 earthquake. 



2. Models for nonuniform fault slip using the RSF law 

The asperity model indicates that spatial heterogene- 
ity of material property is important, and it is com- 
patible with the RSF law discussed in II. C. Regions of 
velocity-weakening frictional property (a — 6 < 0) corre- 
spond to asperities, where stick-slip occurs, and aseismic 
sliding occurs at regions of velocity-strengthening fric- 
tional property (a — 6 > 0). Afterslip occurs in velocity- 
strengthening areas, and it slowly relaxes stress increases 
generated by nearby earthquak es. Using a siiigle-deg ree- 
of-freedom spring-block model, Im arone aLl (|l991l ) ob- 
tained theoretical slip time function u(t) of afterslip, 
which occurs on a fault with velocity-strengthening fric- 
tion (a — 6 > 0), as follows: 



(g - b)a„ 
u(t) = ; m 



kVr> 



(a — b)a, 



-t+l 



+ Vbt, (58) 



where (t„ is normal stress on the fault plane, k is spring 
stiffness, Vcs is coseismic slip velocity, Vq is preseismic 
slip rate, time t is measured from the earthquake occur- 
rence time. Quantitative comparison between afterslip 
observations and models indicate that th e RSF law well 
expla ins afterslip (jFreedl . 12Q07; Pcrfctti ni and Avouad . 
1200411 . 

In case the stiffness is larger than the critical stiffness 
defined by Eq. for a velocity-weakening fault, it 

is called conditionally stable (Scholz, 1988). Although 
aseismic sliding usually occurs under quasi-static loading 
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FIG. 25: Spatial distribution of interplate coupling 
estimated from geodetic data (colored circles) along the 
Sunda trench, where the Australian plate subducts 
beneath the Sumatra island. Red and orange circles 
indicate that the plate interface is nearly locked and 
strain is accumulated during an interseismic period, and 
white and yellow circles indicate that continuous 
aseismic sliding occurs and strain is not accumulated. 
Red and green contours with 5m interval show seismic 
slip in the 2004 Sumatra- Andaman (M = 9.1) and the 
2005 Nias-Simeulue (M = 8.7) earthquakes. Blue and 
black lines show the approximate source areas of the 
1797 and 1833 great earthquakes (|Chlieh et aLl . l2008l) . 



for conditionally stabl e case, rapid st ress increase may 
generate seismic slip ( Gu et all . Il984| ). This fact indi- 
cates that sliding behavior at a fault is not determined 
only by the fault properties but by a loading condition, 
suggestive of variable sliding behavior of a fault. Note 
that the effective stiffness of a fault may be related to 
fault size as will be shown in the next subsection. 

Since the RSF law takes into consideration time- 
dependent healing pr ocess, it can be used in simulations 
of earthquake cycles. iTse and Ricd (198^ first published 
an earthquake cycle model for a strike-slip fault in an 
elastic continuum using the RSF law to successfully ex- 
plain stick-slip behavior at a shallower part of a fault, 
continuous stable sliding at a deeper part, and after- 
slip at intermediate depths. In the simulation, quasi- 
dynamic equilibrium between frictional stress and elastic 
stress generated by fault slip and relative plate motion is 
numerically solved. Their assumption on depth depen- 
dence of a — & is consistent with laboratory data, which 
indicate a—b changes from negative to positive at about 



300° C (iBlanpied et aLl . ll995D . Similar models have been 
presented for earthquake cycles at particular regions to 
compare the simulations with observed earthquake recur- 
rence and/or crustal deformation. Figure UHl shows an 
example simulation result of spatiotemporal evolution of 
slip velocity on a model plate interface, where great in- 
terplate earthquakes repeatedly occ ur at a shallower part 
and stable sliding on a deeper part (jllori e t at, 2004). 

If a single asperity exists on a fault plane without any 
interactions with other asperities, regular stick-slip at a 
constant recurrence interval is expected to occur. Note 
that when the asperity size is close to the critical nu- 
cleation zone size, irregular stick- slip cycle is o bserved 
even for a single asperity model ( Liu and Ricd . [2 007') . 
When some asperities exist within short distances, they 
interact with each other, resulting in complex earthquake 
sequences including single asperity ruptures and multi- 
ple asperity ruptures. Numerical simulations of complex 
earthquake sequences due to interac tions between some 
asperit ies have b een c ar ried out by Kato and Hirasawal 



(Il999bl). iKatol (|2004|) . iLapusta and Liul (l2009l ). and 



iKaneko et al. M2OIOD . In these studies, friction obeying 
the RSF law was assumed and different values of friction 
parameters (a',&',£) are assigned for model asperities 
with velocity- weakening friction, to reproduce compound 
earthquakes, where some asperities are ruptured simulta- 
neously or wit h som e time delays, which resembles some 
observations. iKatd (I2008D . for instance, reproduced a 
complex earthquake cycle similar to that observed at the 
Sanriku-oki region, northeastern Japan, where simulated 
earthquakes included the 1968 Tokachi-oki earthquake 
(M=8.2), the 1994 Sanriku-oki earthquake (M=7.7) and 
its largest aftershock (M=6.9) and afterslip. These stud- 
ies suggest that spatial distribution of asperities or fric- 
tion parameters controls regularity and complexity of 
earthquake recurrence. This further suggests that numer- 
ical forecasts of earthquakes may be possible if we can ob- 
tain detailed map of friction parameters on a fault. Fric- 
tion parameters have actually been estimated through 
compar ison of observ e d dat a and simul ations at Cali- 
forni a (iJohnson et all 20061) and Japan ( Fukuda et all . 
120091: iMivazaki et all l200l from afterslip data. 

Preseismic sliding, which is aseismic sliding during a 
slip nucleation process, is expected before earthquake 
occurrence from the RSF law. It is almost ubiqui- 
tously observed in laboratory experiments, where the 
amo unt of preseismic sliding is of the order of microme- 
ters ( Ohnaka and Shenl . Il999l ). Using a spring- block sys- 
tem implemented with the RSF law, one can show that 
th e preseismic sliding amount is approximately given by 
C (iPopov et a/.l . l2010l ). Some model studies with the RSF 
law discussed crustal deformation expec ted from pre- 
seism i c sliding for parti c ular e arthquakes ( Kuroki et all . 
I2OO2I: IStuart and Tulhd . Il995[) . However, it is difficult 
to predict precise amplitudes of crustal deformation, be- 
cause friction parameters that influence preseismic slid- 
ing are not well constrained from presently available data. 
There are some reports of observations of preseismic 
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FIG. 26: Snapshots of simulated slip rate V on the model plate interface normalized by the relative plate velocity 
Vpi in a model for recurrence of great earthquakes along the Nankai trough, central Japan. Red, white, and blue 
show seismic slip rates, stable sliding with sliding velocity n early equal to the plate velocity, and nearly locked, 

respectively. Modified from lHori et al\ ( 20041 ). 



sliding, thoug h insignifican t or questionable observations 
are included (|Wvssl . ll997t ). For example, the close and 
dense geodetic observation of the Parkfield segment of 
the San Andreas fault could not detect any precursory 
slip prior to the 2004 earthquake, although it should be 
remarked that an observation of the tremor may sug- 
gest the accelerated creep on the fault ^ 16 km be - 
neath the eventual ea r thqua ke hypoccnter ( Shellvl . [2009( ). 
iKanamori and Cipail ( 1974h detected precursory signals 
in long-period strain seismogram before the occurrence 
of the 1960 great Chilian earthquake (M=9.5). Since no 
earthquake that could explain the observed strain sig- 
nals was detected, they inferred that the signals were 
caused by preseismic slid ing on a deeper extensio n of the 
mainshock fault plane. iLinde and Sacksl ( 20021) exam- 
ined crustal deformation data before the occurrence of 
the 1944 Tonankai (M=8.0) and 1946 Nankai (M=8.1) 
earthquakes, southwestern Japan to construct a model 
for preseismic sliding of these earthquakes. Their model 
indicates that preseismic sliding took place at a deeper 
extension of the main shock fault plane. However, in 
models with the common RSF law, accelerating preseis- 
mic sliding just before earthquake occurrence takes place 
within the source area of seismic slip because sponta- 
neous accelerating slip can be nucleated only in velocity- 
weakening region, be ing inconsist ent with these models 
of preseismic sliding. iKatol ( 2003al ) proposed a model for 
earthquake cycles at a subduction zone to explain large 
preseismic sliding at the deeper extension of the scis- 
mogenic plate interface. He assumed velocity-weakening 
friction {diJ,ss/d\nV < 0) at low velocities and velocity- 
strengthening friction {dfiss/d\nV > 0) at high veloci- 
ties, where /Xss is a steady-state friction coefficient given 
by Eq. ((T9| . Preseismic sliding relaxes regional stresses, 
which may decrease seismic activity, while it increases 
stresses around the edges of the sli pped region , whic h 
tends to increase seismic activity ( Kato et all Il997t) . 



This may explain precursor y seismic quiescence observed 
for so me large earthquakes ( Kanamorl Il98ll : IWvss et al\ . 
Il98l[ ). Preseismic shding perturbs regional stress field, 
resulting in an increase or decr ease of seismicity. Tak- 
ing into consideration this effect lQgatal (I2005D systemat- 
ical searched seismicity changes in Japan to find possible 
crustal stress changes due to preseismic sliding. 



3. Slow earthquakes 

Slow earthquakes are episodic fault slip events that gen- 
erate little or no seismic waves because their source du- 
rations are longer than the periods of observable seis- 
mic waves. Slip events without seismic wave radiation 
are often called silent earthquakes or slow slip events. 
Slow earthquakes have been stu died by using records of 
very- long-period seismographs ( Kanamori and Steward . 
Il979[ ). creepmeters that direct l y det ect fault creep at 
the ground s u rface ( King et fflZ.I . [l973l ). and strainmeters 
(|Linde et all Il996l ). Afterslip and preseismic sliding 
mentioned earlier may be included in slow earthquakes. 

Recent development of dense geodetic observation net- 
works including GPS and bore hole tiltmeters acceler- 
ates s tu dies of slow earthqua kes ( Schwartz and Rokoskvl 
l2007t ). iHirose et all ( 1999t ) found that episodic aseis- 
mic slip with duration of about 300 days took place in 
1997 on the plate boundary in the Hyuganada region, 
southwestern Japan, from GPS data. The estimated 
slip and source area indicated that it released seismic 
moment corresponding to magnitude of 6.6. Later, al- 
most the same size aseismic slip events occurred at the 
same area in 2003 and 2010. In the Tokai region, cen- 
tral Japan, another large slow earthquake from 2000 to 
2005 released seismi c moment nearl y equa l to that of an 
M=7 .0 earthquake ( Mivazaki et all 120061 : lOzawa et all 
I2OO2I) . The source area of this slow earthquake was esti- 
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mated at the deeper extension of the locked plate bound- 
ary, where a magnitude 8 class interplate earthquake is 
expected to occur. At almost the same area, smaller 
slow earthquakes, corresponding to moment magnitude 
of about 6.0, with shorter du rations of a few days wer e 
found to occur repeatedly ( Hirose and Obarai l2006( ). 
These slow earthquakes are often called short-term slow 
slip events (SSEs) to be discriminated from long-term 
SSEs of durations of several mon ths or longer. Fur- 
thermore, [Hirose^^nT^Obar^ ( 20061 ) found low-frequency 
tremors, which radiate seismic waves with long dura- 
tions, from high-sensitivity borehole seismometer array 
data. These events are clearly distinguished from long 
durations of wave trains and lack of high-frequency com- 
ponents of seismic waves. Short-term SSEs and low- 
frequency tremors occur simultaneously almost at the 
same locations. Synchronized occurrence of short-term 
SSEs and low-frequency tremors were observed in other 
regions s uch as the Cascadia subd uction zone. North 
America ( Rogers and Drager^ 20031) and Shikoku, south- 
western Japan ( Obara et aLl . [2004 ). 



These findings of slow earthquakes and low-frequency 
tremors force us to reconsider simple view of earthquakes 
as brittle fracture. Many mechanical models for slow 
earthquakes has been proposed. Since both seismic and 
aseismic slip can be easily modeled with the RSF law, 
it is natural to consider that slow earthquakes can be 
modeled with the RSF law. In fact, sustaining aseismic 
oscillation, similar to recurrence of slow earthquakes, oc- 
curs in a single-degree-of-freedom spring-block model if 
the spring stiffness k is equal to the critical stiffness k^^t 
given by Eq. ([22]) (iRuinal Il983l) . 



^ ^ ^ ^ ^ .. _ _ _ Usi ng a more realistic 

elastic continuum model. iKatol ( 2003 ) showed that slow 
earthquakes occur when the size of velocity-weakening 
region is close to the critical size of slip nucleation zone. 
An effective stiffness kes of a fault may be defined by 



fcoff = At/Au, 



(59) 



wher e At is shear-st ress change on the fault due to slip 
Au (iDieterichl . Il986l ). For a circular fault of radius r 
with a constant stress drop in an infinite uniform elastic 
medium with Poission ratio = 0.25, fccfi is given by 



IttG 

krfi = , 

24 r ' 



(60) 



where G is rigidity. Recalling that unstable slip occurs 
for k < fccrt for a spring-block model, unstable slip is 
expected to occur for fcoff < fccrt on a fault in an elastic 
medium. This leads to the condition for occurrence of 
unstable slip is that the fault radius r is larger than the 
crtical fault size rc given by 



G 



24 [b' ~ a')a, 



(61) 



where dn is the normal stress. Note that the critical 
nucleation zone size r^, obtained by considering the sta- 
bility around steady-state sliding may not be realistic 



in natural conditions during earthquake cycles. Other 
forms of critical nucleation zone sizes were obtai n ed by 
considering more re a listic conditions (|Dieterichl . Il992l: 
iRubin and Ampuerol . l2005f ). It is confirmed in numer- 
ical simulations that usual earthquakes with short slip 
duration occurs for a circular fault with r > Tc, contin- 
uous stable sliding for r ^ Tc, and slow earthquakes for 
r ^ rc, where slip duration increases with a decrease in 
r/rc as shown in Fig [27l (Kato, 20q3b, 2 004). The same 
idea was adopted bv lLiu and Ricd (l200iD in their model 
for slow earthquakes at a subduction zone, where they 
showed that high pore fluid pressure in the fault zone 
is required to explain observed recurrence intervals and 
slip amounts of slow earthquakes. Although these mod- 
els are simple and plausible, slow earthquakes may occur 
under limited conditions of r r^. This seems to be 
inconsistent with the observations that slow earthquakes 
are common phenomena at some reg ions. Us ing a two- 
degree -of- freedom spring-block model, lYoshida a nd Kat3 
(I2003D showed that slow earthquakes may occur for wider 
conditions by considering interaction between unstable 
block where usual earthquakes repeatedly occur and a 
conditionally stab l e bloc k wh ere slow earthquakes occur. 
IShibazaki and liol ( 20031 ) and IShibazaki and Shimamotol 
( 20071 ) introduced a cut-off velocity to the state evolution 
effect in Eq. (fTO|) to obtain frictional property of veloc- 
ity weakening [dpLss/dhiV < 0) at low velocities and of 
velocity strengthening {di^iss/dlnV > 0) at high ve loci- 
tics, which is similar to the model bv iKatol (20033) for 
deep preseismic sliding. Similar complex frictional be- 
havior that d^ss/dlnV depends on ve locity was a c tually 
observed in the lab oratory for halite |Shimamotal . Il986l ) 
and for serpentine (jMoore et a/.l . ll997l ). In this case, slip 
is accelerated at low velocities with dfiss/dlnV < and 
is decelerated at high velocities with dfiss/dlnV > 0, 
leading to slow earthquakes. Repeating slow earthquakes 
at transition depths from shallow lock ed zone to deeper 
stable slidi ng zone were simulated in IShibazaki and liol 
(l2003ll andihi bazaki and Shimamotol (|2007D . This kind 
of model was further extended to simulate short- and 
long-term SSEs and the ir interaction with shallo wer large 
interplate earthquakes ( Matsuzawa et aLl . l2010[ ). Weak- 
ness of these models is that experimental data for velocity 
dependence of dfj,ss/dlnV are insufficient and frictional 
properties at depths w here slow earthquakes occur are 
unknown. iRubinI I 2008f ) reviewed models for slow earth- 
quakes based on the RSF law and pointed out that the 
existing models seem to be difficult to explain common 
occurrence of slow earthquakes at subduction zone s. He 
suggested that the variation of pore fluid pressure due 
to inelastic dilation of fault zone and fluid diffusion is 
required for generating slow earthquakes. 
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FIG. 27: The duration of slip events versus r/vc in 
numerical simulations using the RSF law, where a 
circular asperity of radius r with velocity-weakening 
friction is embedded in a fault of velocity-strengthening 
frictional property. Vc denotes the crit ical fault rad ius 
for unstable slip defined in the text ( Katol . [2003bf ) . 



4. Origin of complexities of earthquakes and aftershock decay 
law 



iRicd (I1993D claimed that complex earthquake sequences 
simulated in inherently discrete models may be arti- 
fact and geometrical and/or material heterogeneity is re- 
quired to explain observed complexity of earthquakes. 
Continuum models with relatively homogeneous fric- 
tional properties produce simple patterns of earthquakes 
such as periodic recurrence of large earthquakes that 
break the entire scismogenic zone. Using a c ontin - 
uum model with the RSF law, iBen-Zion and Ric^ ( 19951 ) 
introduced heterogeneity in effective normal stress on 
the fault and successfully produces moderately complex 
earthquake sequences. They pointed out that abrupt 
change in the effective norma l stress is nec e ssary to pro- 
duce complex earthquakes. iHillers et all ( 2007t) intro- 
duced spatial heterogeneity in characteristic slip distance 
C in the model of a vertical strike-slip fault to produce 
complex earthquake sequences that include a wide range 
of earthquake magnitude. The obtained relation be- 
tween earthquake magnitude and frequency mimics the 
Gutenberg-Richter (GR) law and the statistical proper- 
ties of simulated earthquakes depend on the degree of 
heterogeneity in C They also argued temporal cluster- 
ing of simulated ea rthquakes and tendency o f nuclcation 
sites of smaller C. iHiUers and MiUeil (|2007l ) introduced 
spatial variation of pore pressure to generate complex 
earthquake sequences. 

An important fact about the relation between magni- 
tude and frequency of earthquakes obtained in observa- 
tions is that the GR law may not always be valid for 
each individual fault. For some faults and plate bound- 
aries, the number of small earthquakes is too few than 



that expected from the GR law and the frequency of 
large earthquakes that rupt ure the entire fault, indicat- 
ing violation of the GR law ( Ishibe and Shimazakil . 120091 : 
IStirling et al. I. I1996D . This behavior of fewer small earth- 
quakes than that expected from the frequency of large 
earthquakes is referred to the characteristic earthquake 
model. Highly coupled plate interface in the Tokai re- 
gion, central Japan, is nearly quiescent, while many small 
earthquakes o ccur in the overridi ng plate and subducting 
oceanic plate ( Matsumural . Il997t) . This suggests that ex- 
cept for great earthquakes few small earthquakes occur 
at the plate interface in the Tokai region. Considering 
large earthquakes along the San Andreas fault, Califor- 
nia, and smaller e arthquake s at se condary faults around 
the San Andreas, iTurcott j ( 19971 ) argued that the ob- 
served GR law comes from a fractal distribution of faults 
and characteristic earthquakes at each fault. 

Another important empirical law that demonstrates 
complexities of earthquakes is the modified Omori 
(Omori-Utsu) law for decay in aftershock occurrence rate 
( Utsu et aUliggsl ). Aftershock rate n at time t from the 
occurrence of the mainshock is well approximated by 



n{t) 



K 



{t + Tmol)^ 



(62) 



where, K, TmoL; and p arc constants. Constant p takes 
1 for many cases. In the case of Tmol = 0, this re- 
lation is simply referred to as the Omori law. After- 
shocks have been thought to be manifestation of relax- 
ation of stress generated by the mainshock. To explain 
delay times of aftershocks, subcritical crackin g due to 
stress corrosion ( Yamashita and Knopofj . Il987() and the 
variation of effective normal stress due to diffusion of 
pore fluid, who s e pres sure is perturbed by the ma i nshock 
(jBosl and Nu3 . |2002| ). were invoked. iDieterichI (Il994l ) 
considered responses of many fault patches, where fric- 
tion is assumed to obey the RSF law, to instantaneous 
stress change due to the mainshock. He furtehr assumed 
that a constant seismicty rate is achieved under a con- 
stant loading rate without any stress perturbation. This 
model successfully explains the power law decay of seis- 
micity rate with p = 1, being consistent with observa- 
tions, and has been applie d to analyses of a ftershocks 
of some large earthquakes ( Toda et al\ . Il998l ). Another 
important model for aftershocks using the RSF law is 
related to afterslip. Afterslip perturbs stresses around 
its source area, causing aftershocks. Differentiating the 
slip function Eg. ([55)) of afterslip with respect to time, we 
have a slip rate approximately proportional to (t -\- c)~^ , 
which m ay related to a stress rate and therefore scismic- 
ity rate i Perfettini and Avouad . |2004[) . This expected 
scisniicity rate coincides with the Omori-Utsu formula 
with p — 1. Moreover, afterslip propagates outward 
from a main shock slip ar ea, leading to expansion of after- 
shock area l|Katol . l2007t) . Aftershock expansion pattern 
obtained from a numerical model with the RSF law is 
consistent with observed expansions of aftershock ares 
(|Peng and Zhacl . 120091: [Taiima and Kanamoril [l985l) . 
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FIG. 28: Schematic diagram of the relation between 
frictional force and shp distance during shp- weakening 
process, where Fi, Fp, Fj_, and denote the initial 
force, the peak frictional force, the dynamic frictional 
force, and the critical slip distance, respectively. The 
shaded area indicates the fracture energy Gc- 



5. Earthquake dynamics: critical slip distance 

Here we consider the dynamics of unstable motion. The 
unstable slip of a spring-block system given by Eq. ((2T|) 
is accompanied by the drop of frictional force. If one 
plots the frictional force as a function of the slip distance 
(Fig. one can define the distance Dc over which 

the frictional force drops. This behavior of decreasing 
frictional force wit h increasing slip i s ref e rred t o the slip- 
weakening model (lAndrewsl . llQTfil: lldal . Il972[ ). and the 



slip distance Dc is called the critical slip distance in seis- 
mology. Dc is on the same order of (or at most several 
tens of) the characteristic length C in an evolution law 
( Bizzarri and Coccd . l2003D . This is so irrespective of the 
number of degrees of freedom: discrete or continuum. 

Importantly, one can estimate Dc of earthquakes by 
analyzing seismic wave. Such analyses show that Dc is 
on the order of sever al tens of centimeters or a meter 
( Ide and Taked . 119971) . Note that fracture energy Gc, 
which is equal to twice the surface energy density F, can 
rather stably be estimated from seismic waveform data, 
though accurate estimate of Dc is difficult because of 
poor re solution of rupture process m odeling from seismic 
waves ( Guatteri and Spudishl . [2OO0I) . The characteristic 
slip distance C estimated for afterslip of a large interplate 
earthquake by GPS data is on the order of mm (Fukuda 
et al. 2009). This makes a quite contrast to laboratory 
experiments, where C is typically estimated as several 
micrometers. Because £ is a typical longitudinal dimen- 
sion of true contact patch, application of the RSF law to 
a natural fault implies that a natural fault also consists 
of true contact patches, a linear dimension of which is 
several tens of centimeters. Although the aperture of a 
fault is not empty but filled with fluid and gouge, a fault 
generally has the non-planer structure (e.g., jogs) that 
can interlock to resist the displacement. Such jogs may 
effectively act as the true contact area. However, it is 
not obvious at all if the RSF law still holds for such true 



contact area of macroscopic scale. 

At least, we believe that the RSF law should not be 
used except for very low speed friction. Namely, the RSF 
law no longer holds at seismic slip rate due to physical 
processes caused by frictional heat: flash heating, melt- 
ing, mechanochcmical effects, etc. In such cases, the crit- 
ical slip distance Dc is proportional to ec/ P, where ec is 
the critical energy per unit area for a weakening process 
(e.g. melting) to occur and P is the normal stress. (As 
the frictional force is proportional to P, the produced 
heat is proportional to P and to the slip distance D. 
Thus, the weakening process may occur if PD is on the 
order of ec.) Namely, the critical slip distance is inversely 
proportional to the normal pressure. This implies that 
the critical slip distance must be smaller for deeper faults. 
However, unfortunately, such depth dependence has not 
been observed so that the mechanism that determines 
the critical slip distance is different. 

Another important process that affect the critical slip 
distance is the off-fault fracture accompanied by the 
crack propagation on fault. Andrews (2005) analyzes a 
model for the slip propagation on a fault supplemented 
with the Coulomb yield condition for off-fault material. 
He flnds that the effective critical slip distance depends 
on the distance from the crack initiation point. This is 
because the plastic zone is wider for larger crack. Thus, 
the critical slip distance is essentially scale dependent, 
which is consistent with the observation facts. 



IV. EARTHQUAKE MODELS AND STATISTICS II: SOC 
AND OTHER MODELS 

A. Statistical properties of the OFC model 

1. The model 

In the previous section, we reviewed the properties of sta- 
tistical physical models of earthquakes such as the spring- 
block BK model and the continuum model. In the present 
section, we deal with further simplified statistical phys- 
ical models of earthquakes (Turcotte 1997; Hergarten, 
2002; Turcotte, 2009). Many of them were coupled map 
lattice models originally introduced as the SOC models 
(Bak, Tang and Wiesenfeld, 1987; Bak and Tang, 1989, 
Ito and Matsuzaki, 1990; Nakanishi, 1990; Brown, Scholz 
and Rundle, 1991; Olami, Feder and Christensen, 1992; 
Hainzl, Zoller and Kurths, 1999; 2000; Hergarten, and 
Neugebauer, 2000; Helmstetter, Hergarten and Sornette, 
2004). 

The one introduced by Olami, Feder and Christensen 
(OFC) as a further simplification of the BK model, now 
called the OFC model, is particularly popular (Olami, 
Feder and Christensen, 1992). It is a two-dimensional 
coupled map lattice model where the rupture propa- 
gates from lattice site to its nearest-neighboring sites 
in a non-conservative manner, often causing multi-site 
"avalanches" . Extensive numerical studies have also been 
devoted to this model, mainly in the field of statistical 
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physics, which we wish to review in the present section 
(Christensen and Olami, 1992; Grassberger, 1994; Mid- 
dleton and Tang, 1995; Bottani and Delamotte, 1997; 
de Carvalho and Prado, 2000; Pinho and Prado, 2000; 
Lise and Paczuski, 2001; Miller and Bouher, 2002; Her- 
garten and Neugebauer, 2002; Boulter and Miller, 2003; 
Helmstettcr, Hergarten and Sonette, 2004; Pcixoto and 
Prado, 2006; Wisscl and Drossel, 2006; Ramos, 2006; 
Kotani, Yoshino and Kawamura, 2008; Kawamura et al, 
2010, Jagla, 2010). 

In the OFC model, "stress" variable fi {fi > 0) is as- 
signed to each site on a square lattice with L x L sites. 
Initially, a random value in the interval [0,1] is assigned 
to each fi, while fi is increased with a constant rate uni- 
formly over the lattice until, at a certain site i, the fi 
value reaches a threshold, /c = 1. Then, the site i "top- 
ples", and a fraction of stress afi (0 < a < 0.25) is 
transmitted to its four nearest neighbors, while fi itself 
is reset to zero. If the stress of some of the neighboring 
sites j exceeds the threshold, i.e., fj > fc~ 1, the site 
j also topples, distributing a fraction of stress afj to its 
four nearest neighbors. Such a sequence of topplings con- 
tinues until the stress of all sites on the lattice becomes 
smaller than the threshold fc- A sequence of toppling 
events, which is assumed to occur instantaneously, cor- 
responds to one seismic event or an avalanche. After an 
avalanche, the system goes into an interseismic period 
where uniform loading of / is resumed, until some of the 
sites reaches the threshold and the next avalanche starts. 

The transmission parameter a measures the extent of 
non-conservation of the model. (This a should not be 
confused with a describing the velocity-weakening fric- 
tion force employed in the study of the BK model of 
subsection III. A. We are using a as a conservation pa- 
rameter of the OFC model throughout this subsection 
IV. A). The system is conservative for a = 0.25, and is 
non-conservative for a < 0.25. A unit of of time is taken 
to be the time required to load / from zero to unity. 

In the OFC model, boundary conditions play a cru- 
cial role. For example, SOC state is realized under open 
or free boundary conditions, but is not realized under 
periodic boundary conditions. Thus, most of the stud- 
ies made in the past employed open (or free) boundary 
conditions. 



2. Properties of the homogeneous model 

Earlier studies concentrated mostly on the event size dis- 
tribution of the model (Olami, Fedcr and Christensen, 
1992; Christensen and Olami, 1994; Grassberger, 1994; 
de Carvalho and Prado, 2000; Lise and Paczuski, 2001; 
Miller and Bouher, 2002; Bouher and Miller, 2003; 
Drossel, 2006). The avalanche size s is defined by the 
total number of "topples" in a given avalanche, which 
could be larger than the number of toppled sites because 
multi-toppling is possible in a given avalanche. (In fact, 
it is observed that multi-toppling rarely occurs in the 



model except in the conservation limit or in the regime 
very close to it.) It turned out that the size distribution 
of the model exhibited a power-law-like behavior close 
to the GR law. Yet, there still remains controversy con- 
cerning whether the model is strictly critical (Lise and 
Paczuski, 2001) or only approximately so (de Carvalho 
and Prado, 2000; Miller and Boulter, 2002; Bouher and 
Miller, 2003; Drossel, 2006). 

In Fig[^ we show the size distribution of the model 
under open boundary conditions for several values of the 
transmission parameter a (Kawamura et al, 2010). As 
can be seen from the figure, a near straight-line behavior 
corresponding a power-law is observed. The slope repre- 
senting the B-value is not universal varying from ~ 0.90 
to ~ 0.22 as a is varied from 0.17 to 0.245. The power-law 
feature is weakened as one approaches the conservation 
limit. 
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FIG. 29: (Color online) The size distribution of the 
OFC model under open boundary conditions for various 

values of the transmission parameter a. The slope of 
the data gives the value of 1 + i?, which is shown in the 
figure. Taken from (Kawamura et al, 2010). 

Hergarten et al observed that the OFC model also 
exhibited another well-known power-law feature of seis- 
micity, i.e., the Omori law (or the inverse Omori law) 
describing the time evolution of the frequency of after- 
shocks (foreshocks) (Hergarten and Neugebauer, 2002; 
Helmstetter, Hergarten and Sonette, 2004). We show in 
FiglSUTa) on a log-log plot the frequency of aftershocks 
as a function of the time elapsed after the mainshock 
t (Kawamura et al, 2010). The slope representing the 
Omori exponent p is again not universal depending on the 
parameter a as p = 0.84, 0.69 and 0.03 for a = 0.17,0.20 
and 0.23, respectively. Since the p- value is known to come 
around unity in real seismicity, the p value of the OFC 
model is not necessarily close to real observation. Simi- 
lar results are obtained also for foreshocks: See FiglSOTb). 
Aftershocks and foreshocks are defined here as events of 
arbitrary sizes which occur in the vicinity of mainshock 
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with its epicenter lying with the distance r < Vc (the 
range parameter it taken to be = 10 in FiglHU]). 
As one approaches the conservation hmit a ~ 0.25, both 
aftershocks and foreshocks tend to go away. 
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FIG. 30: (Color online) The time dependence of the 
frequency of aftershocks (a), and of foreshocks (b), of 
the OFC model under open boundary conditions on a 

log-log plot for several values of the transmission 
parameter a. Mainshocks are the events of their size 
greater than s > ^ 100. The time t is measured with 
the occurrence of a mainshock as an origin. The range 
parameter is Vc — 10. Taken from (Kawamura et al, 
2010). 

As mentioned, the properties of the model depends 
on applied boundary conditions. Middleton and Tang 
observed that the model under open boundary condi- 
tions went into a special transient state where events 
of size 1 (single-site events) repeated periodically with 
period l — 4a (Middleton and Tang, 1995). These single- 



site events occur in turn in a spatially random manner, 
but after time 1 — 4a, the same site topples repeatedly. 
Although such a periodic state consisting of single-site 
events is a steady state under periodic boundary con- 
ditions, it is not a steady state under open boundary 
conditions because of the boundary. Indeed, clusters are 
formed near the boundary, within which the stress val- 
ues are more or less uniform, and gradually invades the 
interior destroying the periodic state. Eventually, such 
clusters span the entire lattice, giving rise to an SOC- 
like steady state. Middleton and Tang pointed out that 
such clusters might be formed via synchronization be- 
tween the interior site and the boundary site, the latter 
having a slower effective loading rate due to the bound- 
ary. Large-scale synchronization occurring in the steady 
state of the OFC model was further investigated by Bot- 
tani and Delamottc (Bottani and Dclamotte, 1997). 



In contrast to the aforementioned critical properties of 
the model, recent studies also unraveled the character- 
istic features of the OFC model (Ramos, 2006; Kotani, 
Yoshino and Kawamura, 2008; Kawamura et al, 2010). 
By investigating the time series of events, Ramos found 
the nearly periodic recurrence of large events (Ramos, 
2006). Kotani et al studied the spatiotemporal correla- 
tions of the model and identified in the OFC model a phe- 
nomenon resembling the "asperity" (Kotani, Yoshino and 
Kawamura, 2008; Kawamura et al, 2010). These authors 
computed the local recurrence-time distribution, P{T), 
of the model. The computed P{T), shown in Fig. [211 ex- 
hibited a sharp (5-function-like peak at T = T* = 1 — Aa, 
indicating that many (though not all) events of the OFC 
model were repeated with a fixed time-interval T = T*. 
While the peak at T = T* is sharp, it is not infinitely 
sharp with a finite intrinsic width: See the inset. The 
peak position turned out to be independent of the range 
parameter Tc, the size threshold Sc, and the lattice size (as 
long as it was not too small). As a is increased toward 
a ~ 0.25, the (5-function peak is gradually suppressed 
with keeping its position strictly at T = 1 — 4q;. The S- 
function peak of P{T) goes away toward the conservation 
limit a = 0.25: See Fig.[3TJ 



In the longer time regime T > T* , P[T) exhibits be- 
haviors close to power laws (Kotani, Yoshino and Kawa- 
mura, 2008; Kawamura et al, 2010). Furthermore, the 
periodic events contributing to a sharp peak of P(T) 
( "peak events" ) possess a power-law- like size distribution 
very much similar to those observed for other aperiodic 
events (Kotani, Yoshino and Kawamura, 2008; Kawa- 
mura et al, 2010). Hence, in earthquake recurrence of 
the model, the characteristic or periodic feature, i.e., a 
sharp peak in P{T) at T = T* , and the critical feature, 
i.e., powcr-law-likc behaviors of P{T) at T > T* and 
power-law- like size distribution, coexist. 



38 




FIG. 31: (Color online) Log-log plots of the local 
recurrence-time distributions of large avalanches of their 
size s > Sc = 100 for a fixed range parameter rc ~ 10, 
with varying the transmission parameter a. The arrow 
in the figure represents the expected peak position for 

a = 0.245 corresponding to the period 
T* = 1 — Aa = 0.02. The inset is a magnified view of 
the main peak for the case of a = 0.17. Taken from 
(Kawamura et al, 2010). 



3. Asperity-like phenomena 

In fact, it has turned out that the (5-function peak of P{T) 
is borne by "asperity-like" events, i.e., the events which 
rupture repeatedly with almost the same period 1 — 4a 
and with a common rupture zone and a common epicen- 
ter. In seismology, the concept of asperity is now quite 
popular. A typical example might be the one observed 
along the subduction zone in northeastern Japan, partic- 
ularly repeating earthquakes off Kamaishi (Matsuzawa, 
Igarashi and Hasegawa, 2002; Okada, Matsuzawa and 
Hasegawa, 2003). 

In Fig. [221 we show typical examples of such asperity- 
like events as observed in the OFC model (Kawamura 
et al, 2010). In the upper panel, we show for the case 
of a = 0.2 typical snapshots of the stress- variable distri- 
bution immediately before and after a large event which 
occurs at time t ^ to. Discontinuous drop of the stress 
associated with a rupture of a synchronized cluster is 
discernible. Then, at time t = Iq + T* , the same cluster 
(except for a minor difference) ruptures again. In the 
lower panel, we show snapshots of the stress- variable dis- 
tribution immediately before and after this subsequent 
avalanche occurring at t ^ Iq +T* . In this particular ex- 
ample, a rhythmic rupture of essentially the same cluster 
has repeated more than ten times. 

The asperity-like events go away in the conservation 
limit a — >■ 1/4 (Kawamura et al, 2010). It is also ob- 
served that an epicenter site tends to lie at the tip or at 
the corner of the rupture zone rather than in its interior 
(Kawamura et al, 2010). The asperity-like events ob- 



served in the OFC model closely resemble those familiar 
in seismology (Scholz, 2002), in the sense that almost the 
same spatial region ruptures repeatedly with a common 
epicenter site and with a common period. 




FIG. 32: (Color online) Snapshots of the stress-variable 
distribution of the OFC model under open boundary 
conditions for the case of a = 0.2; (a) immediately 
before a large event at time t = to, (b) immediately 
after this event, (c) immediately before the following 
event which occurs at time t = to + T*{T* — 0.2), and 
(d) immediately after this second event. Two events are 
of size s = 15891 and s = 15910 on a i = 256 lattice. 
The region surrounded by red bold lines represents the 

rupture zone, while the star symbol represents an 
epicenter site which is located at the tip of the rupture 
zone. Taken from (Kawamura et al, 2010). 

In fact, not all large events of the OFC model oc- 
cur in the form of asperity. Many clusters forming large 
events are left out of the rhythmic recurrence, and rup- 
ture more critically with widely-distributed recurrence 
time, thereby bearing the observed power-law-like part 

of p(r). 

A key ingredient in the asperity formation is a self- 
organization of the highly concentrated stress state 
(Kawamura et al, 2010). The stress- variable distribution 
in the asperity region tends to be "discretized" to cer- 
tain values. In Fig. [331 we show for the case of a = 0.17 
the stress- variable distribution D{f) of the asperity sites 
immediately (a) before and (b) after an avalanche, aver- 
aged over asperity events. As can be seen from the figure, 
D{f) now consists of several "spikes" located at appro- 
priate multiples of the transmission parameter a, i.e., 
at 1 — na before the rupture, and at / = na after the 
rupture, with n being an integer. Furthermore, as the as- 
perity events repeat, the tendency of the stress-variable 
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concentration is more and more enhanced. In Fig. l34l 
we show the time sequence of the stress-variable distri- 
bution at the time of topphng for the asperity events. 
As the asperity events repeat, the stress- variable distri- 
bution tends to be narrower, being more concentrated 
on the threshold value /c = 1 (Kawamura et al, 2010; 
Hcrgartcn and Krenn, 2011). 

In fact, one can prove that the stress- variable distri- 
bution at the time of toppling tends to be more con- 
centrated on the threshold value /c = 1 as the asper- 
ity events repeat (Kawamura et al, 2010). Namely, once 
each site starts to topple more or less at similar stress 
values close to the threshold value /c = 1, this ten- 
dency is more and more evolved as the asperity events 
repeat. The stress-variable concentration tends to he self- 
organized. Such a stress- variable concentration imme- 
diately explains why the interval time of the asperity 
events is equal to 1 — Aa, and why the same site becomes 
an epicenter in the asperity sequence (Kawamura et al, 
2010). For example, the reason why the interval time 
is 1 — 4a when all sites topple at the stress value close 
to the threshold /c = 1 in the asperity events can easily 
be seen just by remembering the conservation law of the 
stress, i.e., the stress- variable dissipated at the time of 
toppling, which is 1 — 4a per site if the toppling occurs 
exactly at / = 1, should match the stress loaded during 
the interval time T. See (Kawamura et al, 2010), for fur- 
ther details. More recently, Hergarten and Krenn made 
further analysis of this stress-concentration phenomenon, 
demonstrating that the mean stress excess representing 
the extent of the stress concentration approaches zero ex- 
ponentially with a certain decay time which is dependent 
on the number of "internal" sites (the sites contained in 
the rupture zone) connected to an epicenter site (Her- 
garten and Krenn, 2011). Then, the epicenter site with 
the smallest number of internal nearest-neighbor sites, 
i.e., the one lying at the tip of the rupture zone, has 
the longest decay time and turns out to be most stable. 
This observation gives an explanation of the finding of 
(Kawamura et al, 2010) that the majority of epicenter 
sites of the asprity-like events are located at the tip of 
the rupture zone. 

Although the origin of the asperity is usually ascribed 
in seismology to possible inhomogeneity of the material 
property of the crust or of the external conditions of 
that particular region, wc stress that, in the present OFC 
model, there is no built-in inhomogeneity in the model 
parameters nor in the external conditions. The "asper- 
ity" in the OFC model is self- generated from the spatially 
uniform evolution-rule and model parameters. 

As mentioned, the asperity in the OFC model is not 
a permanent one: In long terms, its position and shape 
change. After all, the model is uniform. Nevertheless, re- 
covery of spatial uniformity often takes a long time, and 
the asperity exists stably over many earthquake recur- 
rences. Although one has to be careful in immediately 
applying the present result for the OFC model to real 
earthquakes, it might be instructive to recognize that the 
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FIG. 33: (Color online) The stress- variable distribution 
D{f) of each site contained in the rupture zone of the 
asperity event of the OFC model under open boundary 

conditions, just before (a) and after (b) the asperity 
event. An asperity event is defined here as an event of 
its size greater than s > Sc = 100 belonging to the main 
peak of the local recurrence-time distribution function. 
The transmission parameter is a = 0.17. The inset is a 
magnified view of the main peak. Taken from 
(Kawamura et al, 2010). 



observation of asperity-likc earthquake recurrence does 
not immediately mean that the asperity region possesses 
different material properties nor different external condi- 
tions from other regions. 

Thus, critical and characteristic features coexist in the 
OFC model in an intriguing manner. Although the criti- 
cal features were emphasized in earlier works, the model 
certainly involves the eminent characteristic features in 
it as well. Thus, the OFC model, though an extremely 
simplified model, may capture some of the essential ingre- 
dients necessary to understand apparent coexistence of 
critical and characteristic properties in real earthquakes. 
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FIG. 34: (Color online) The time sequence of the 
stress- variable distribution D{f) at the time of toppling 

of each site contained in the rupture zone of the 
asperity events. An asperity event is defined here as an 
event of its size greater than s > Sc = 100 belonging to 
the main peak of the local recurrence-time distribution 
function. The transmission parameter is a = 0.17. As 
the events repeat, the stress- variable distribution at the 
time of toppling gets more and more concentrated on 
the borderline value /c = 1. Taken from (Kawamura et 
al, 2010). 



4. Effects of inhomogeneity 

It should be noticed that the original OFC model is a 
spatially homogeneous model, where homogeneity of an 
earthquake fault is implicitly assumed. Needles to say, 
real earthquake fault is more or less spatially inhomoge- 
neous, which might play an important role in real seis- 
micity. Then, a natural next step is to extend the origi- 
nal homogeneous OFC model to the inhomogeneous one 
where the evolution rule and/or the model parameters 
arc taken to be random from site to site. 

Spatial inhomogeneity could be either static or dynam- 
ical. As a cause of possible temporal variation of spatial 
inhomogeneity, one may consider the two distinct pro- 
cesses, i.e., the fast dynamical process during an earth- 
quake rupture changing the fault state via, e.g., wear, 
frictional heating, melting, etc and many slower processes 
taking place during a long interseismic period until the 
next earthquake, e.g., water migration, plastic deforma- 
tion, chemical reactions, etc (Scholz, 2002). Thus, in in- 
troducing the spatial inhomogeneity into the OFC model, 
there might be two extreme ways: In one, one may as- 
sume that the randomness is quenched in time, namely, 
spatial inhomogeneity is fixed over many earthquake re- 
currences. In the other extreme, spatial inhomogeneity 
is assumed to vary with time in an uncorrelated way over 
earthquake recurrences. 

Several studies have been made on the inhomogeneous 
OFC model for both types of inhomogencities. For the 
first type of inhomogeneity, i.e., the quenched or static 



randomness, Janosi and Kertesz introduced spatial in- 
homogeneity into the stress threshold and found that 
the inhomogeneity destroyed the SOC feature of the 
model (Janosi and Kertesz, 1993). Torvund and Froyland 
studied the effect of spatial inhomogeneity in the stress 
threshold, and observed that the inhomogeneity induced 
a periodic repetition of system-size avalanches (Torvund 
and Froyland, 1995). Ceva introduced defects associated 
with the transmission parameter a, and observed that 
the SOC feature was robust against small number of de- 
fects (Ceva, 1995). Mousseau and Bach et al introduced 
inhomogeneity into the transmission parameter at each 
site. These authors observed that the bulk sites fully syn- 
chronized in the form of a system-wide avalanche over a 
wide parameter range of the model (Mousseau , 1996; 
Bach, Wisscl and Dresscl, 2008). 

For the second type of inhomogeneity, i.e., the dy- 
namical randomness, Ramos considered the randomness 
associated with the stress threshold, and observed that 
the nearly periodic recurrence of large events persisted 
(Ramos, 2006). More recently, Jagla studied the same 
stress-threshold inhomogeneity, to find that the GR law 
was weakened by randomness (Jagla, 2010). Very inter- 
esting observation by Jagla is that, once the slow struc- 
tural relaxation process is added to the inhomogeneous 
OFC model, both the GR law and the Omori law are 
realized with the exponents which are stable against the 
choice of the model parameter values and are close to the 
observed values. Yamamoto et al studied the dynami- 
cally inhomogeneous model with a variety of implemen- 
tations of the form of inhomogeneities to find the gen- 
eral tendency that critical features found in the original 
homogeneous OFC model, e.g., the Gutenberg- Richter 
law and the Omori law, were weakened or suppressed 
in the presence of inhomogeneity, whereas the character- 
istic features of the original homogeneous OFC model, 
e.g., the near-periodic recurrence of large events and the 
aspcrity-likc phenomena, tended to persist (Yamamoto, 
Yoshino and Kawamura, 2010). 

Thus, the properties of the dynamically inhomoge- 
neous models are quite different from those of the static 
or quenched inhomogeneous models. In the latter case, 
introduced inhomogeneity often gives rise to a full syn- 
chronization and a periodic repetition of system-size 
events. Such a system-wide synchronization is never 
realized in the dynamically homogeneous models. Pre- 
sumably, temporal variation of the spatial inhomogene- 
ity may eventually average out the inhomogeneity over 
many earthquake recurrences, giving rise to the behavior 
similar to that of the homogeneous model. 



B. Fiber bundle models 



The fiber bundle model, initiated bv lPeircd(|l926D in the 
context of testing strength of cotton yarns, represents 
various aspects of fracture processes of disordered sys- 
tems, through its self-organised dynamics (for detailed 
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FIG. 35: The fiber bundle consists initially of A'' fibers 
attached in parallel to a fixed and rigid plate at the top 
and a downwardly movable platform from which a load 
W is suspended at the bottom. In the equal load 
sharing model considered here, the platform is 
absolutely rigid and the load W is consequently shared 
equally by all the intact fibers. 



review see lPradhan et al. 1 (|2010D '). The fiber bundle (see 
Fig. [55]) consists of N fibers or Hook springs, each having 
identical spring constant k. The bundle supports a load 
W = Na and the breaking threshold {crth)i of the fibers 
are assumed to be different for different fiber {i). For 
the equal load sharing model we consider here, the lower 
platform is absolutely rigid, and therefore no local de- 
formation and hence no stress concentration occurs any- 
where around the failed fibers. This ensures equal load 
sharing, i.e., the intact fibers share the applied load W 
equally and the load per fiber increases as more and more 
fibers fail. The strength of each of the fiber {(7th) i in the 
bundle is given by the stress value it can bear, and be- 
yond which it fails. The strength of the fibers are taken 
from a randomly distributed normalised density p{(Jth) 
within the interval and 1 such that 



p{(Tth)dath = 1- 



The equal load sharing assumption neglects 'local' fluc- 
tuations in stress (and its redistribution) and renders the 
model as a mean-field one. 

The breaking dynamics starts when an initial stress 
a (load per fiber) is applied on the bundle. The fibers 
having strength less than a fail instantly. Due to this 
rupture, total number of intact fibers decreases and rest 
of the (intact) fibers have to bear the applied load on 
the bundle. Hence effective stress on the fibers increases 
and this compels some more fibers to break. These two 
sequential operations, namely the stress redistribution 
and further breaking of fibers continue till an equilib- 
rium is reached, where either the surviving fibers are 



FIG. 36: The simple model considered here assumes 
uniform density p{crth) of the fiber strength distribution 

up to a cutoff strength (normalized to unity). At any 
load per fiber level ut at time t, the fraction at fails and 
1 — CTt survives. 



strong enough to bear the applied load on the bundle 
or all fibers fail. 

This breaking dynamics can be represented by recur- 
sion relations in discrete time steps. For this, let us con- 
sider a very simple model of fiber bundles where the fibers 
(having the same spring constant k) have a white or uni- 
form strength distribution p{uth) upto a cutoff strength 
normalized to unity, as shown in Fig. 1361 p{(Jth) = 1 for 
< <7th < 1 and p{(Jth) = for ath > 1- Let us also de- 
fine [/((cr) to be the fraction of fibers in the bundle that 
survive after (discrete) time step t, counted from the time 
t = when the load is put (time step indicates the num- 
ber of stress re-distributions). As such, Ut{a = 0) = 1 for 
all t and Ut{cr) = 1 for f = for any a; Ut{a) = U*{a) ^ 
for i — > oo and a < af, the critical or failure strength of 
the bundle, and Ut{<j) = for i — >■ oo if cr > ct/. 

Therefore Ut{(j) follows a simple recursion relation (see 
Fig. 



t+i 



I -at] at = 



W 



or, Ut+i = 1 



a 



At the equilibrium state {Ut+i = Ut 
relation takes a quadratic form of U* 



(63) 

C/*), the above 



U* 



U* 



The solution is 



U*{a) = -±{af-a) 



cr = 0. 



Here aj \s the critical value of initial applied stress be- 
yond which the bundle fails completely. The solution 
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with (+) sign is the stable o ne, whereas the one wi t h (— ) 
sign gives unst a ble so l ution (iBhattacharvva et al\ . 20031 : 
iPradhan etai\ . l2002t iPradhan and ChakrabartiL l200lh . 
The quantity U* (cr) must be real valued as it has a phys- 
ical meaning: it is the fraction of the original bundle that 
remains intact under a fixed applied stress a when the 
applied stress lies in the range < cr < cr/. Clearly, 
t/*(0) = 1. Therefore the stable solution can be written 



U*{a) = U*{af) + {af - af'^- U*{af) = \ and = \. 

, (64) 

For 17 > <Tf wc can not get a real- valued fixed point as 
the dynamics never stops until Ut = Q when the bundle 
breaks completely. 



At the critical point (cr = cr/), we observe a different 
dynamic critical behavior in the relaxation of the failure 
process. From the recursion relation it can be shown 
that decay of the fraction Ut {af) of unbroken fibers that 
remain intact at time t follows a simple power-law decay: 



Ut 



t + i 



(69) 



starting from Uq ^ I. For large t (t — > oo), this reduces 
to Ut — 1/2 oc (5 = 1; a strict power law which is a 
robust characterizat ion of the critical state (see, however, 
IZapperi et ail (|l997t )). 



1. Universality class of the model 



(a) At a < af 

It may be noted that the quantity U*{a) — U*[(jf) be- 
haves like an order parameter that determines a transi- 
tion from a state of partial failure (cr < cr/) to a state of 
total failure (cr > cr/): 



OEEC/*(a)-C/*(a/) = (a/-a)'3;/?: 



(65) 



To study the dynamics away from criticality {a ^ af 
from below), wc replace the recursion relation (|63p by a 
differential equation 



dU 
'~dt 



-U + (7 

u ' 



Close to the fixed point we write C/t(cr) — U*{a) + e 
(where e 0). This, following Eq. (p5)) . gives 

e = C/t((T) - J7*(cr) « cxp(-t/T), (66) 

l] . Near the critical point 



whereT= i [i(a/ -a)-i/2 
we can write 



(67) 



Therefore the relaxation time diverges following a power- 
law as CT —> cr/ from below. 

One can also consider the breakdown susceptibility x, 
defined as the change of U*{cr) due to an infinitesimal 
increment of the applied stress a 



X 



dU*{a) 



da 



i(f7/-a)-^;7=i (68) 



from Eq. EH Hence the susceptibility diverges as the 
applied stress cr approaches the critical value cr/ = i. 
Such a divergence in x had already been observed in the 
numerical studies. 



(b) At a = af 



The universality class of the model has been checked tak- 
ing two other types of fiber strength distributions: (I) 
linearly increasing density distribution and (II) linearly 
decreasing density distribution within the {uth) limit 
and 1. One can show that while cr/ changes with differ- 
ent strength distributions (cr/ = ■\/4/27 for case (I) and 
CT/ ~ 4/27 for case II), the critical behavior remains un- 
changed: a = 1 /2 = = 7, ^ = 1 for all these equal load 
sharing models ( Pradhan et a/.l . [20101 ). 



2. Precursors of global failure in the model 

In any such failure case, it is important to know when 
the failure will take place. In this model, there exist 
several precursors. The growth of susceptibility x with 
CT, following Eg. indeed suggests one such possibility: 
Y~^/ ^ decreases linearly with cr to at cr = cr/ from be- 
low. IPradhan and Hemmerl ( 20091 ) studied the rate R{t) 



(= — ^jf-) of failure of fibers following the dynamics like 
in Eq. [521 for cr > cr / and found that the rate becomes 
minimum at a time t^, which is half of the failure time 
t/ of the bundle. This relation is shown to be indepen- 
dent of the breaking strength di stribution of the fibers . 
A sim ilar relation was also found ( Pradhan and Hemmeil . 

[mil) 

for the rate of energy released in a bundle. This is, 
of course, easier to measure using accoustic emmisions. 



3. Strength of the local load sharing fiber bundles 

So far, we studied models with fibers sharing the exter- 
nal load equally. This type of model shows (both ana- 
lytically and numerically) existence of a critical strength 
(non zero cr/) of the macroscopic bundle beyond which it 
collapses. The other extreme model, i.e., the local load 
sharing model has been proved to be difficult to tackle 
analytically. 

It is clear, however, that the extreme statistics comes 
into play for such local load sharing models, for which 
the strength cr/ —> as the bundle size {N) approaches 
infinity. Essentially, for any finite load (cr), depending on 



43 




FIG. 37: The breaking rate R{t) vs. the rescaled step 
variable tf /t for the uniform threshold distribution for a 
bundle of = 10^ fibers. Different symbols are for 
different excess stress levels a — af : 0.001 (circles), 
0.003 (triangle s), 0.005 (squares) and 0.007 (crosses). 
From ( Pradhan and Hemmerl . |2009| ) . 



the fiber strength distribution, the size of the defect clus- 
ter can be estimated using Lifshitz argument (see section 
III.A[) as InA^, giving the failure strength aj ^ l/(lnA^)°, 
where the expo nent a assumes a value appropriate for the 
model fsee e.g.. iPradhan and Chakrabartil ( 2003bl) ). If a 
fraction / of the load of the failed fiber goes for global re- 
distribution and the rest (fraction 1 — /) goes to the fibers 
neighboring to the failed one, then we see that there is 
a crossover from extreme to self-aver aging statistics at a 
finite value of / (see e.g., |Pradhan~t at, (,2010l) ). 



4. Burst distribution: crossover behavior 

In fiber bundle models, when the load is slowly increased 
until a new failure occurs, a burst can be defined as the 
number (A) of fiber failures following that failure. The 
distribution of such bursts (D(A)) shows power-law be- 
havior. It was shown for a generic case (independent of 
threshold distribution) that the form of this distribution 
(for continuous loading) is 



D{A)/N = CA 



(70) 



for 



in the limit N ^ oo. 

The burst exponent C has a value | 
age over all a(= to af) and it is universal 
( Hemmer and Hansenl . [19921 ) . However, the burst expo- 
nent value depends, for e.g., on the details of loading 
process and also from which point of the loading the 
burst statistics are recorded. If the burst distribution 
is recorded only near the critical point (a S, erf), the 
burst exponent (C) value becomes 3/2 ([Pradhan et all . 
|2005[ ). For equal load sharing model model with uniform 
strength distribution, the burst distribution is shown 
fFig. [55|) for recording that starts from different points of 




FIG. 38: The burst size distribution for different values 
of xq in the equal load sharing model with uniform 
threshold distribu tion. The number of fi bers is 
TV = 50000 (jPradhan et al\ . I2006D . 
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FIG. 39: Crossover signature in the local magnitude 
distribution of earthquakes in Japan. During the 100 
days before mainshock the expone nt is 0.60; much 
smaller than the average value 0.88. ( KawamuraL l2006h 



effective loading which is denoted by xq, where Xt = cr/Ut 
is the elongation or the effective loading (for linear elas- 
tic behavior) at any point t). The crossover behavior is 
clearly seen. In these studies, the load increase rate is 
extremely slow and the increase is assumed to stop once 
a fiber fails. The consequent avalanches are studied at 
that load. Once the avalanche stops, the load is increased 
again. This process is realistic in the case of earthquakes 
where stress accumulation takes place over years. How- 
ever, if the increase in load is fixed (da), then the above 
exponent value of C becomes 3: A ~ '^^^da ^ ' S^^^^S 
A^2 ^ ^ _ (from Eq. [64)) and since D{A)dA ~ da, 
D{A) ~ ~ A-«, C = 3 (IPradhan et aLl . l2002D . 

In fact, the earthquake frequency statistics may indeed 
show the crossover behavior mentioned above: If event 
frequency is denoted by D{M), then it is known that 
D{M) ^ M~^, where M denotes the magnitude (may be 
assumed to be related to avalanche s i ze A in the mod- 
els) and C value is found ( Kawamural l2006l ) to be more 
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FIG. 40: (a) Schematic representation of the rough 
earth's surface and the tectonic plate, (b) The 
one-dimensional projection of the surfaces form 
overlapping Cantor sets. 




FIG. 41: (a) Two Cantor sets along the axes r and 
r — r' . (b) The overlap si(r) along the diagonal, (c) 
The corresponding density pi{s). 



(C « 0.9) for statistics over a smaller time period (before 
the mainshock), compared to the long time average value 
(C w 0.6); sec Fig.|39l 



C. Two fractal overlap model 

The common geometrical property observed in seismic 
faults is its fractal nature. It is now well known that, 
like other fractured surfaces, fault surfaces also p osses 
self-afHne roughness (see e.g.. ISantucci et all (2003) and 
references therein) . Therefore, it is worth investigating if 
earthquake phenomena can be modelled as an outcome 
of relat ive movement of two self-afhne s urfac es over each 
other. IChakrabarti and Stinchcombd ( 1999I ). in a sim- 
plistic model, studied the overlap statistics of two Can- 
tor sets in order to understand the underlying physics of 
such phenomena. 

Cantor set is a prototype example of fractal. In order 
to construct a triadic Cantor set, in the first step the 
middle third of a base interval [0,1] is removed. In the 
successive steps, the middle thirds of the remaining inter- 
vals ([0,1/3] and [2/3,1] and so on) are removed. After n 
such steps, the remaining set is called a Cantor set of gen- 
eration n. When this process is continued ad infinitum 
i.e., in the limit n — >■ oo, it becomes a true fractal. 

In this model, the solid-solid contact surfaces of both 
the earth's crust and the tectonic plate are considered 
as average self-affine surfaces (see Fig. |40]). The strain 
energy grown between the two surfaces due to a stick pe- 
riod is taken to be proportional to the overlap between 
them. During a slip event, this energy is released. Con- 
sidering that such slips occur at intervals proportional 
to the length corresponding to that area, a power-law 
for the frequency distribution of the energy release is ob- 
tained. This compares well with the GR law (see e.g. 
iBhattacharvva and Chakrabartil ( 20061 )). 



1. Renormalisation group study: continuum limit 

Let the sequence of generators Gn define the Cantor set 
at the n-th generation within the interval [0,1]: Gq ~ 
[0, 1], Gi = RGo = [0, a] U [&, 1], ... ,G„+i = i?G„, ... . 
The mass density of the set Gn is represented by -D„(r) 




FIG. 42: The overlap densities (probabilities) p{s) at 
various generations; (a) zeroth, (b) first, (c) second, and 
(d) infinite generation. 



i.e., Dn{r) = 1 if ?' is in any of the occupied intervals of 
Gn and elsewhere. The overlap magnitude between the 
sets at any generation n is then given by the convolution 
form Sn{r) = J dr' Dn(j'')Dn{r — r') (for symmetric frac- 
tals). One can express the overlap integral si in the first 
generation by the projection of the shaded region along 
the vertical diagonals in Fig. l4lT a). That gives the form 
shown in Fig. I41f b). For a = b < 1/3, the non- vanishing 
si(r) regions do not overlap and are symmetric on both 
sides with slope of the middle curve being exactly double 
those on the sides. One can then easily check that the 
distribution pi(s) of overlap s at this generation is given 
by Fig. HH with both c and d greater than unity, main- 
taining the normalisation condition with cd = 5/3. The 
successive generations of the density Pn{s) may therefore 
be represented by Fig. [551 where 



Pn+l(s) = Rpn(s) = -Pn(-) + ^P«( — )• (71) 

5 c 5 c 

In the infinite generation limit of the renormalisation 
group (RG) equation, if p* (s) denotes the fixed point 
distribution such that p*{s) = Rp*{s), then assuming 
p*{s) - s'-'pis), one gets {d/5)c' + {U/b){c/2)^ = 
1. Here p{s) represents an arbitrary modular function, 
which also includes a logarithmic correction for large s. 
This agrees with the above mentioned normalisation con- 
dition cd = 5/3 for the choice 7 = 1 giving 



p*{s) = p{s) 



7 



1 



(72) 



The above analysis is for the continuous relative mo- 
tion of the overlapping fractals. For discrete steps. 
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the contact area distribution can be found exactly for 
two Cantor sets hav ing same dimension (log2/log3) 
( Bhattacharvval 120051 ) . The step size is taken as the min- 
imum element in the generation at which the distribution 
is found. 



2. Discrete limit 

Let Sn{t) represent the amount of overlap between the 
two Cantor sets of generation n at time t. Initially {t = 0) 
the two identical Cantor sets arc placed on top of each 
other, generating the maximum overlap (2" for the n- 
th generation sets). Then in every time step (discrete) 
the length of the shift is chosen to be 1/3" for the n- 
th generation, such that a line segment in one set ei- 
ther completely overlaps with one such segment on the 
other set or does not overlap at all, i.e., partial overlap 
of two segments of the two sets are not allowed. Periodic 
boundary conditions are assigned in both the sets. The 
magnitude of overlap (s„(t)), therefore, in this discrete 
version, is given by the number of overlapping pairs if line 
segment of the two sets. Because of the structure of the 
Cantor sets, the overlap magnitudes can only have cer- 
tain discrete values which are in geometric progression: 
Sn = 2"-^ k = 0,...,n. 

Let Nr(sn) denote the the number of times an overlap 
Sn has occurred in one period of the time series for the 
n-th generation (i.e. 3" time steps). It ca n be shown 
that ( Bhattacharvva and Chakrabartii |2006| ) 



Nr{2"-'') = "Cfc2'' 



fc = 0, . 



(73) 



Now, if Prob{sn) denotes the probability that after time 
t there arc s„ overlapping segments, then for the general 
case of Sn = 2""*^, fe = 0, . . . , n it is given by 

_ iVr(2""^-) 



Prob{2 



J2 Afr(2"-'^) 

fe=0 

2^ 

_ /l\"-^/2^^- 

' ^""Msi I 3, 



(74) 



3. Gutenberg-Richter law 

Since the allowed values of the overlap are s„ = 2"^*^, 
k = 0, . . . , n, one can write log2 Sn = n ~ k. Then the 
above equation becomes 



Probisn) = "Clo 



iog,.„(3y V3 

= F(l0g2 Sn). 

Near the maxima it may be written as 



Tl-log2 S„ 



(75) 



F{M) 



3 . 9(M-n/3)^ 

-^], (76) 



where M = log2 Sn- To obtain the GR law analog from 
this distribution we have to integrate F{M) from M to 
oo. 

Fcura{M) = / F{M')dM' 



M 
oo 



3 , 9(A/'-n/3)^ 

exp(--^ L^)dM'{77) 

Z\/mr 4 n 



Substituting p ~ j^iM' — n/3) we get 

oo 

FcumiM) = J exp{-p^)dp. 

On simplification, it gives 



(78) 



^(A/-«/3) 



1 hi 

Fcum{M) = -W - CXp 
3 V TT 



9 [M-n/'df 



{M - n/i)- 



(79) 

Fcum{M) in the above equation suggests that the 'aver- 
age' quakes arc of magnitude ?t./3, while 

Fcu,n{M) ^ cxp [-(9/4)(Af - n/if/n] (80) 



can be simplified for large M. Using e " 

+ 00 ^ ^ +00 

— oo — oo 

where xo refers to the extremal point with d f / dx\x=xa 

(9/4)[J\/(mo/Vi)--2Af/3] 



l/^/2^) / e-="'/2-i-V2ax^2; and / erf^'^'^dx - e"-'''^'') 



0, one finds Fcum,{M) 
g-3A//4 ^Yiere xq = (j^) ™o; 
(jBhattacharva et all [20091) 



mo 



log F,um{M) = A - -M, 



n. It gives 



(81) 



2\/mr 



where A is a constant depending on n. This is the 
Gutenberg-Richter law in the model and clearly holds 
for the high magnitude end of the distribution. Also, one 
can equate easily the magnitude M with the released en- 
ergy E by noting that M = log2 s here. The overlap s is 
related to energy E and hence the relation M ~ log£^, 
giving Fcum ~ E-^/^. 

Similar to outcome of the simple fractal models con- 
sidered here, a power law behavior for the overlap distri- 
bution also occurs for two overlapping random Cantor 
sets, Sierpinsky gasket and Sier pinsky carpet o verlap- 
ping on their respective replica ( Pradhan et all . [2003i') . 
and a fractional Brownian p rofile overlapping on an- 
other ( De Rubeis et all Il996[ ) . In view of the generality 
of the power law distribution and the fractal geometry 
of the fault surfaces, it is suggested that the GR law 
owes its origin significantly to the fractal geometry of 
the fault surfaces. It may be noted that identifying the 
aftershocks as these adjusted overlaps, with average size 
given by Eq. (j79p . one can define an average magnitude 
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(n/3) dependent on the fractal geometry generator frac- 
tion (= 1/3 here) and the genration number (n). This 
agrees with the ob s erved data quite satisfactorily (see 
iBhattacharva et al\ (|201lD V 



4. Omori law 

Let N^^^°''{t) denote the cumulative number of after- 
shocks (of magnitude M > Mq, where Mq is some thresh- 
old) after the mainshock. Then the Omori law states that 



dt ~ tP' 



(82) 



The value of the exponent p is close to unity for tec- 
tonically active region, alth ough a range of p values a re 
also observed (for review see IBhattacharva et al\ ( 20091 )). 
In practice, a particular value of p is observed when the 
threshold Mq is given. For this model, when the thresh- 
old is fixed at the minimum (i.e., Mq = 1), then p = Q 
due to the fact that aftershock occurs at every step in 
this model. However, interesting facts are seen when the 
threshold is set at the second highest possible value n—1 
(recall that the second highest overlap was 2"~^). Then 
for t = 2.3^^ (where, ri — 0, . . . ,ti — 1) there is an af- 
tershock of magnitude n—1. Therefore, neglecting the 
prefactor 2, an aftershock of magnitude n—1 occurs in 
geometric progression with common ratio 3. Therefore 
we get the general rule N{3t) = N{t) + 1, leading to 



N{t)=\og,it). 



(83) 



On integration, Omori law gives N{t) — t^~P, and there- 
fore from this model we get p = 1, which is the Omori 
law suggested value for p. The model therefore gives a 
range of p values between and 1 which systematically 
increases within the range of threshold values. 



V. DISCUSSIONS AND CONCLUSIONS 

Earthquakes, due to their devastating consequences, have 
been a subjected of extensive studies in various diciplines, 
ranging from seismology to physics. Although the ef- 
feorts were not always commensurate (see also iKaganI 
(l2006l) for a critical view of the inherent difficulties of 
the present approach of theretical physics), in the last 
decade considerable progress have been made in study- 
ing different aspects of this vast topic. In this review, the 
progresses in such studies are discussed from the point of 
view of statistical physics. 

Principally being a large scale dynamic failure process, 
it is necessary to formulate the background of friction and 
fracture in order to understand the physics of earthquake. 
In Sec. II such issues are discussed: After the Griffith's 
theory for crack nuclcation and the fracture stress statis- 
tics of disordered solids, we discuss the RSF law and 



microscopic models for solid-solid friction. Also, the ef- 
fects that could lead to violations of RSF laws are also 
discussed. 

Several statistical approaches to model earthquake dy- 
namics are discussed. The BK model is discussed in one 
and two spatial dimensions as well as its long range ver- 
sion. In Sec. IIIA6, the continuum limit is also discussed, 
which gives 'characteristic' earthquakes. BK model has 
also been discussed in terms of RSF law. Apart from 
relatively complex modelling like that of BK models and 
continuum models, we also discuss simplistic models such 
as OFC models, fiber bundle models and purely geomet- 
rical models like the Two Fractal Overlap model. While 
many details are lost in any such model, they still cap- 
tures the complex nature of the dynamics and the differ- 
ent statistical aspects, helping us to gain new insights. 

As one can easily see that inspite of considerable 
progress in the study of such an important and complex 
dynamical phenomenon as earthquake, our knowledge is 
far short of any satisfactory level. We believe, major col- 
laborative efforts, involving physicist and seismologicsts 
in particular, are urgently necessary to unfold the dynam- 
ics and employ our knowledge of the precursor events to 
save us from catastrophic disasters in future. 
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Appendix A: Glossary 

aftershocks: Small earthquakes that follow a large 
earthquake (main shock). 

afterslip: Aseismic sliding that follows an earthquake, 
asperity: (a) A region where stick-slip motion occurs on 
a fault or a plate boundary. Strain energy is accumulated 
at an asperity during a stick stage between earthquakes 
and it is released by seismic slip at the occurrence of 
an earthquake, (b) Junction of protrusions of the two 
contacting surfaces. 

Cantor set: One starts with the set of real numbers in 
the interval [0 : 1] and divide the set in a few subsets and 
remove one of the subsets in the first step. As the removal 
scheme is repeated ad infinitum, one is left with a dust of 
real numbers called the Cantor set. It is a fractal object, 
characteristic earthquakes: Earthquakes that repeat- 
edly rupture approximately the same segment of a fault. 
The magnitudes and slip distributions of characteristic 
earthquakes are similar to one another, 
dynamical critical phenomena: Critical behaviors, 
which are associated with the dynamical properties of 
the system, rather than the equilibrium properties (e.g., 



47 



thermal transition is Ising model) are called dynamical 
critical phenomena (e.g., depinning transition of a frac- 
ture front, time dependent field induced transitions in 
Isng model etc.). 

fiber bundle model: Originating from texttile engi- 
neering, fiber bundle model is often used as a prototype 
model for fracture dynamics. In its simplest form, it con- 
sists of a large number of fibers or Hookc-springs. The 
bundle hangs from a rigid ceiling and supports, via a 
platform at the bottom, a load. Each fiber has got iden- 
tical spring constants but the breaking stress for each 
differs. Depending on the breaking stress of the fibers, 
the fibers fail and successive failure occur due to load 
redistribution, showing complex failure dynamics, 
fractals: A fractal is a geometrical object having self- 
similarity in its internal structure. 

fractional Brownian profile : Fractional Brownian 
motion (fBm) is a continuous time random walk with 
zero mean. However, the directions of the subsequent 
steps of an fBm are correlated (positively or negatively). 
A fractional Brownian profile is the trajectory of such a 
walk. It is self-similar. 

Gutenberg-Richter (GR) law: The power law 
describing the magnitude-frequency relation of earth- 
quakes. The frequency of earthquakes of its energy (seis- 
mic moment) E decays with E according to oc = 
E~^^~^3^'> where B and b = |i? are appropriate expo- 
nents. 

Hamiltonian: It is essentially the total energy of a sys- 
tem. For a closed system, it would be the sum of kinetic 
and potential energies. 

Omori law: The power law describing the decay of the 
number (frequency) of aftershocks with the time elapsed 
after the mainshock. 

power-law distribution: (Also called 'scale free distri- 
bution') A distribution of the generic form P{x) ^ x" . 
Note that there is no length-scale associated with this 
type of distribution, since a transformation like x x/b 
would keep the functional form unchanged. Obscrvablcs 
(e.g., magnetisation, susceptibility etc.) show power-law 
behavior near criticality. Therefore it is often considered 
as a signature of critical behavior. 

slow earthquakes: Fault slip events that radiate little 
or no seismic wave radiations. Rupture propagation ve- 
locities and slip velocities of slow earthquakes are much 
smaller than those of ordinary earthquakes. Slow earth- 
quakes without seismic wave radiations are often called 
silent earthquakes. 

quenched randomness: The randomness in the sys- 
tem which is not in thermal equilibrium with the same 
reservoir as the system and does not fluctuate are called 
quenched randomness. 

rate-and-state friction (RSF) law: An empirical con- 
stitutive law describing the dynamic friction coefficient 
either at steady states or transient states, 
self-organized criticality (SOC): When the dynamics 
of a system leads it to a state of criticality (where scale 
invariance in time and space are observed) without any 



need of external tuning parameter, the system is said to 
have self-organized to a critical state. This phenomenon, 
where a critical point is an attractor of the dynamics, is 
called self-organized criticality. 

self-similarity and self-affinity: Self similarity refers 
to the property of an object that it is similar (exactly or 
approximately) to one or more of its own part(s). Self- 
affinity refers to the properties of those objects which, 
in order to be self-similar, are to be scaled by different 
factor in x and y direction (for 2-d object). 
Sierpinski gasket and Sierpinski carpet: Sicrpinski 
carpet is a fractal object, embedded in a 2-d surface. Its 
construction is as follows: First a square is taken and it 
is divided into 9 equal squares. Then the square in the 
middle is removed, then similar operation is performed 
upon the 8 remaining squares. This process is continued 
ad infinitum to obtain what is called a Sierpinski carpet. 
Sierpinski gasket (also called Sierpinski triangle) is again 
a fractal object. Its construction is as follows: First a 
equilateral triangle is taken. Then it is divided into four 
equilateral triangle of same sizes and the middle one is 
removed. Then same operation is performed upon the 
three remaining triangles. When this process is continued 
ad infinitum, one is left with what is called the Sicrpinski 
gasket. 

universality class: Phase transitions are characterised 
by a set of critical exponent values. The values of these 
exponents are independent of the microscopic details of 
the system and only depend on the symmetry and di- 
mensionality of the order parameter. Therefore, a large 
class of systems often have same critical exponent values. 
A Universality class is a group of systems having same 
critical exponent values. 
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